From a1986d5655a85df377a40cdd53a724efaf1c61b6 Mon Sep 17 00:00:00 2001 From: taw27 Date: Tue, 19 Feb 2008 21:20:25 +0000 Subject: Initial simplex refinement stuff git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@262 bf6ca9ba-c028-0410-8290-897cf20841d1 --- src/refine.c | 417 ++++++++++++----------------------------------------------- 1 file changed, 86 insertions(+), 331 deletions(-) (limited to 'src/refine.c') diff --git a/src/refine.c b/src/refine.c index b18e5fb..45aec8e 100644 --- a/src/refine.c +++ b/src/refine.c @@ -19,394 +19,149 @@ #include #include #include -#include -#include -#include +#include "control.h" #include "displaywindow.h" -#include "gtk-valuegraph.h" -#include "basis.h" -#include "reflections.h" #include "image.h" #include "reproject.h" -#include "control.h" #include "mapping.h" -#include "imagedisplay.h" -#include "utils.h" -typedef enum { - AXIS_X = 1, - AXIS_Y = 2 -} Axis; +/* A simplex is an array of ten of these */ +typedef struct { + double dax; double dbx; double dcx; + double day; double dby; double dcy; + double daz; double dbz; double dcz; +} SimplexVertex; + +typedef struct { + signed int h; signed int k; signed int l; + double dx; double dy; double dz; +} Deviation; + +void refine_do_sequence(ControlContext *ctx) { + + +} -static void refine_fix_unconstrained(gsl_matrix *M) { +void refine_do_cell(ControlContext *ctx) { - int row, j; - int modified = 0; + SimplexVertex s[10]; + Deviation *d; + const double delta = 0.1e9; + int i, nf, f, v_worst; + double fom_worst; - /* Find a row which is all zero */ - for ( row=0; rowsize1; row++ ) { + if ( !ctx->cell_lattice ) { + displaywindow_error("No reciprocal unit cell has been found.", ctx->dw); + return; + } - int zero = 1; + if ( ctx->images->n_images == 0 ) { + displaywindow_error("There are no images to refine against.", ctx->dw); + return; + } + + /* Initialise the simplex */ + s[0].dax = 0.0; s[0].dbx = 0.0; s[0].dcx = 0.0; + s[0].day = 0.0; s[0].dby = 0.0; s[0].dcy = 0.0; + s[0].daz = 0.0; s[0].dbz = 0.0; s[0].dcz = 0.0; + memcpy(&s[1], &s[0], sizeof(SimplexVertex)); s[1].dax = delta; + memcpy(&s[2], &s[0], sizeof(SimplexVertex)); s[2].dbx = delta; + memcpy(&s[3], &s[0], sizeof(SimplexVertex)); s[3].dcx = delta; + memcpy(&s[4], &s[0], sizeof(SimplexVertex)); s[4].day = delta; + memcpy(&s[5], &s[0], sizeof(SimplexVertex)); s[5].dby = delta; + memcpy(&s[6], &s[0], sizeof(SimplexVertex)); s[6].dcy = delta; + memcpy(&s[7], &s[0], sizeof(SimplexVertex)); s[7].daz = delta; + memcpy(&s[8], &s[0], sizeof(SimplexVertex)); s[8].dbz = delta; + memcpy(&s[9], &s[0], sizeof(SimplexVertex)); s[9].dcz = delta; + + /* Create the table of indicies and deviations */ + nf = 0; + for ( i=0; iimages->n_images; i++ ) { + int j; - for ( j=0; jsize2; j++ ) { - if ( gsl_matrix_get(M, row, j) != 0.0 ) { - zero = 0; - break; - } + if ( !ctx->images->images[i].rflist ) { + ctx->images->images[i].rflist = reproject_get_reflections(&ctx->images->images[i], ctx->cell_lattice); } - if ( zero ) { - - /* Find a column which is all zero */ - int i, col; - int row_altered = 0; - - for ( col=0; colsize2; col++ ) { - - int zero2 = 1; - - if ( row_altered ) break; /* This row is now fixed, so move on to the next */ - - for ( i=0; isize1; i++ ) { - if ( gsl_matrix_get(M, i, col) != 0.0 ) { - zero2 = 0; - break; - } - } - - if ( zero2 ) { - gsl_matrix_set(M, row, col, 1.0); - modified = 1; - row_altered = 1; - } - - } - + for ( j=0; jimages->images[i].rflist->n_features; j++ ) { + if ( ctx->images->images[i].rflist->features[j].partner != NULL ) nf++; } } + printf("RF: There are %i partnered features in total\n", nf); - if ( modified) printf("RF: M was modified to fix unconstrained variables.\n"); - -} - -/* Use the IPR algorithm to make "cell" fit the given image */ -ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) { - - ImageFeatureList *flist; - int j, ns; - Axis axis; - double dax = 0.0, dbx = 0.0, dcx = 0.0; - double day = 0.0, dby = 0.0, dcy = 0.0; - double dlx = 0.0, dly = 0.0, dlz = 0.0; - - flist = image->features; - - for ( axis = AXIS_X; axis <= AXIS_Y; axis++ ) { - - gsl_permutation *perm; - int s; - double det; - gsl_matrix *M; - gsl_vector *q; - gsl_vector *p; - gsl_vector *scratch; - gsl_matrix *V; - gsl_vector *S; + d = malloc(nf*sizeof(Deviation)); + f = 0; + for ( i=0; iimages->n_images; i++ ) { - if ( axis == AXIS_X ) printf("RF: ------------------------------------------------------------------ Refining x coordinates -----\n"); - if ( axis == AXIS_Y ) printf("RF: ------------------------------------------------------------------ Refining y coordinates -----\n"); + ImageRecord *image; + int j; - M = gsl_matrix_alloc(3, 3); - p = gsl_vector_alloc(3); - gsl_matrix_set_zero(M); - gsl_vector_set_zero(p); + image = &ctx->images->images[i]; - ns = 0; - for ( j=0; jrflist->n_features; j++ ) { + for ( j=0; jimages->images[i].features->n_features; j++ ) { - double val; ImageFeature *rf; - signed int h, k, l; - double xy; double dix, diy, dx, dy; + double dlx, dly, dlz; double old_x, old_y; rf = &image->rflist->features[j]; - if ( !rf->partner ) continue; - h = rf->reflection->h; - k = rf->reflection->k; - l = rf->reflection->l; + d[f].h = rf->reflection->h; + d[f].k = rf->reflection->k; + d[f].l = rf->reflection->l; /* Determine the difference vector */ dix = rf->partner->x - rf->x; diy = rf->partner->y - rf->y; - printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, h, k, l, dix, diy); + printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, d[f].h, d[f].k, d[f].l, dix, diy); old_x = rf->partner->x; old_y = rf->partner->y; rf->partner->x = dix + rf->partner->parent->x_centre; rf->partner->y = diy + rf->partner->parent->y_centre; mapping_scale(rf->partner, &dx, &dy); + mapping_rotate(dx, dy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); rf->partner->x = old_x; rf->partner->y = old_y; - printf("=> %+10.5f %+10.5f nm^-1\n", dx/1e9, dy/1e9); - - xy = 0; - switch ( axis ) { - case AXIS_X : xy = dx; break; - case AXIS_Y : xy = dy; break; - } + printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/1e9, dly/1e9, dlz/1e9); - /* Elements of "p" */ - val = gsl_vector_get(p, 0); val += xy * h; gsl_vector_set(p, 0, val); - val = gsl_vector_get(p, 1); val += xy * k; gsl_vector_set(p, 1, val); - val = gsl_vector_get(p, 2); val += xy * l; gsl_vector_set(p, 2, val); - gsl_matrix_get(M, 2, 2); - /* Elements of "M" */ - val = gsl_matrix_get(M, 0, 0); val += h * h; gsl_matrix_set(M, 0, 0, val); - val = gsl_matrix_get(M, 0, 1); val += k * h; gsl_matrix_set(M, 0, 1, val); - val = gsl_matrix_get(M, 0, 2); val += l * h; gsl_matrix_set(M, 0, 2, val); + d[f].dx = dlx; + d[f].dy = dly; + d[f].dz = dlz; - val = gsl_matrix_get(M, 1, 0); val += h * k; gsl_matrix_set(M, 1, 0, val); - val = gsl_matrix_get(M, 1, 1); val += k * k; gsl_matrix_set(M, 1, 1, val); - val = gsl_matrix_get(M, 1, 2); val += l * k; gsl_matrix_set(M, 1, 2, val); - - val = gsl_matrix_get(M, 2, 0); val += h * l; gsl_matrix_set(M, 2, 0, val); - val = gsl_matrix_get(M, 2, 1); val += k * l; gsl_matrix_set(M, 2, 1, val); - val = gsl_matrix_get(M, 2, 2); val += l * l; gsl_matrix_set(M, 2, 2, val); - - ns++; + f++; } - if ( ns == 0 ) { - printf("RF: No partners found\n"); - gsl_matrix_free(M); - gsl_vector_free(p); - return NULL; - } - - /* Prepare for solving */ - refine_fix_unconstrained(M); - // gsl_matrix_set(M, 2, 0, 0.0); - // gsl_matrix_set(M, 2, 1, 0.0); - // gsl_matrix_set(M, 2, 2, 1.0); - // gsl_matrix_set(M, 0, 2, 0.0); - // gsl_matrix_set(M, 1, 2, 0.0); - // gsl_vector_set(p, 2, 0.0); - matrix_vector_show(M, p, "RF: "); - - /* Calculate and display the determinant */ - perm = gsl_permutation_alloc(M->size1); - gsl_linalg_LU_decomp(M, perm, &s); - det = gsl_linalg_LU_det(M, s); - printf("RF: Determinant of M = %f\n", det); - gsl_permutation_free(perm); - - /* Solve the equations */ - V = gsl_matrix_alloc(M->size2, M->size2); - S = gsl_vector_alloc(M->size2); - scratch = gsl_vector_alloc(M->size2); - gsl_linalg_SV_decomp(M, V, S, scratch); - q = gsl_vector_alloc(3); /* This is the "answer" */ - gsl_vector_set_zero(q); - gsl_linalg_SV_solve(M, V, S, p, q); - gsl_vector_free(scratch); - gsl_matrix_free(V); - gsl_vector_free(S); - gsl_matrix_free(M); - gsl_vector_free(p); - - switch ( axis ) { - case AXIS_X : { - dax = gsl_vector_get(q, 0); - dbx = gsl_vector_get(q, 1); /* These are the deviations, in the direction "x" of the image coordinate */ - dcx = gsl_vector_get(q, 2); /* system, of a,b and c */ - break; - } - case AXIS_Y : { - day = gsl_vector_get(q, 0); - dby = gsl_vector_get(q, 1); /* These are the deviations, in the direction "y" of the image coordinate */ - dcy = gsl_vector_get(q, 2); /* system, of a,b and c */ - break; - } - } - - gsl_vector_free(q); - } + assert( f == nf ); - printf("RF: ------------------------------------------------------------------ Refinement results ---------\n"); - printf("RF: a should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dax/1e9, day/1e9); - printf("RF: b should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dbx/1e9, dby/1e9); - printf("RF: c should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dcx/1e9, dcy/1e9); - - /* Update the cell */ - mapping_rotate(dax, day, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); - printf("RF: a changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9); - cell->a.x += dlx; cell->a.y += dly; cell->a.z += dlz; - - mapping_rotate(dbx, dby, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); - printf("RF: b changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9); - cell->b.x += dlx; cell->b.y += dly; cell->b.z += dlz; + /* Find the least favourable vertex of the simplex */ + v_worst = 0; + fom_worst = 0; + for ( i=0; i<10; i++ ) { - mapping_rotate(dcx, dcy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); - printf("RF: c changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9); - cell->c.x += dlx; cell->c.y += dly; cell->c.z += dlz; - - return NULL; - -} - -static int refine_sequence_sweep(ControlContext *ctx, double *fit, double *warp) { - - int i; - double series_dev_max = 0; - double series_dev_min = +HUGE_VAL; - double series_dev_mean = 0; - int series_dev_n = 0; - - /* Ensure that ctx->cell_lattice is set up */ - reproject_cell_to_lattice(ctx); - - for ( i=0; iimages->n_images; i++ ) { - - ImageRecord *image; - int j, n; - double image_dev_mean = 0; - - image = &ctx->images->images[i]; + double fom = 0; + int j; - /* Fit this image and update ctx->cell_lattice, index the selected pattern */ - if ( !image->rflist ) image->rflist = reproject_get_reflections(image, ctx->cell_lattice); - refine_fit_image(ctx->cell, image); - reproject_cell_to_lattice(ctx); - image->rflist = reproject_get_reflections(image, ctx->cell_lattice); + for ( j=0; jrflist->n_features; j++ ) { + double xdev, ydev, zdev; - double dix, diy; + xdev = 0; + ydev = 0; + zdev = 0; - /* Skip if no partner */ - if ( !image->rflist->features[j].partner ) continue; + fom += sqrt(xdev*xdev + ydev*ydev + zdev*zdev); - /* Determine the difference vector */ - dix = image->rflist->features[j].partner->x - image->rflist->features[j].x; - diy = image->rflist->features[j].partner->y - image->rflist->features[j].y; - - image_dev_mean += sqrt(dix*dix + diy*diy); - n++; - - } - image_dev_mean /= n; - - if ( image_dev_mean > series_dev_max ) series_dev_max = image_dev_mean; - if ( image_dev_mean < series_dev_min ) series_dev_min = image_dev_mean; - series_dev_mean += image_dev_mean; - series_dev_n++; - - } - - series_dev_mean /= series_dev_n; - *fit = series_dev_mean; - *warp = (series_dev_max - series_dev_min)/series_dev_min; - - return 0; - -} - -void refine_do_stack(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; - - fit_vals = malloc(401*sizeof(double)); - warp_vals = malloc(401*sizeof(double)); - idx = 0; - - if ( !ctx->cell_lattice ) { - displaywindow_error("No reciprocal unit cell has been found.", ctx->dw); - return; - } - - if ( ctx->images->n_images == 0 ) { - displaywindow_error("There are no images to refine against.", ctx->dw); - return; - } - - 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_sweep(ctx, &fit, &warp) ) { - printf("RF: Sequencer sweep failed\n"); - } - 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)); } - displaywindow_update(ctx->dw); - 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); - -} - -void refine_do_image(ControlContext *ctx) { - - if ( !ctx->cell_lattice ) { - displaywindow_error("No reciprocal unit cell has been found.", ctx->dw); - return; - } - - if ( ctx->images->n_images == 0 ) { - displaywindow_error("There are no images to refine against.", ctx->dw); - return; - } - - ImageFeature *fitted; - - fitted = refine_fit_image(ctx->cell, &ctx->images->images[ctx->dw->cur_image]); - ctx->images->images[ctx->dw->cur_image].rflist = NULL; reproject_lattice_changed(ctx); displaywindow_update(ctx->dw); -- cgit v1.2.3