diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2008-02-19 21:20:25 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2008-02-19 21:20:25 +0000 |
commit | a1986d5655a85df377a40cdd53a724efaf1c61b6 (patch) | |
tree | a69caacf764e0272367e29d70666c8f18af06d69 /src | |
parent | af605c813075afa76157f871710116451bfc357a (diff) |
Initial simplex refinement stuff
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@262 bf6ca9ba-c028-0410-8290-897cf20841d1
Diffstat (limited to 'src')
-rw-r--r-- | src/Makefile.am | 5 | ||||
-rw-r--r-- | src/displaywindow.c | 12 | ||||
-rw-r--r-- | src/image.c | 2 | ||||
-rw-r--r-- | src/image.h | 2 | ||||
-rw-r--r-- | src/iprtest.c | 180 | ||||
-rw-r--r-- | src/refine.c | 417 | ||||
-rw-r--r-- | src/refine.h | 9 |
7 files changed, 98 insertions, 529 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index d51f2bb..e9aae1c 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,13 +1,10 @@ -bin_PROGRAMS = dtr iprtest +bin_PROGRAMS = dtr dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c cache.c \ itrans-threshold.c itrans-zaefferer.c itrans-stat.c control.c mapping.c reproject.c prealign.c basis.c \ dirax.c image.c refine.c gtk-valuegraph.c intensities.c glbits.c dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas -lutil -iprtest_SOURCES = iprtest.c reflections.c basis.c utils.c reproject.c image.c refine.c mapping.c -iprtest_LDADD = @LIBS@ @GTK_LIBS@ -lm -lgsl -lgslcblas -lutil - AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@ AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" diff --git a/src/displaywindow.c b/src/displaywindow.c index af17934..5bd00fc 100644 --- a/src/displaywindow.c +++ b/src/displaywindow.c @@ -339,12 +339,12 @@ static gint displaywindow_image_last(GtkWidget *widget, DisplayWindow *dw) { } static gint displaywindow_refinestep(GtkWidget *widget, DisplayWindow *dw) { - refine_do_image(dw->ctx); + refine_do_cell(dw->ctx); return 0; } static gint displaywindow_refinestack(GtkWidget *widget, DisplayWindow *dw) { - refine_do_stack(dw->ctx); + refine_do_sequence(dw->ctx); return 0; } @@ -392,8 +392,8 @@ static void displaywindow_addmenubar(DisplayWindow *dw) { { "DirAxAction", "dtr-dirax", "Start _DirAx", "<Ctrl>D", NULL, G_CALLBACK(displaywindow_dirax) }, { "DirAxReRunAction", NULL, "Run another DirAx cycle", NULL, NULL, G_CALLBACK(displaywindow_dirax_rerun) }, { "StopDirAxAction", NULL, "Stop DirAx", NULL, NULL, G_CALLBACK(displaywindow_dirax_stop) }, - { "RefineStepAction", "dtr-refine", "Refine Current Image", NULL, NULL, G_CALLBACK(displaywindow_refinestep) }, - { "RefineSeqAction", "dtr-refineall", "Refine Entire Stack", NULL, NULL, G_CALLBACK(displaywindow_refinestack) }, + { "RefineStepAction", "dtr-refine", "Refine Unit Cell", NULL, NULL, G_CALLBACK(displaywindow_refinestep) }, + { "RefineSeqAction", "dtr-refineall", "Run Refinement Sequence", NULL, NULL, G_CALLBACK(displaywindow_refinestack) }, { "SetAxisAction", "dtr-tiltaxis", "Set Tilt Axis Position...", NULL, NULL, G_CALLBACK(displaywindow_setaxis) }, { "IncrAxisAction", NULL, "Increase Tilt Axis Position", "<Ctrl>Up", NULL, G_CALLBACK(displaywindow_incraxis) }, { "DecrAxisAction", NULL, "Decrease Tilt Axis Position", "<Ctrl>Down", NULL, G_CALLBACK(displaywindow_decraxis) }, @@ -404,8 +404,8 @@ static void displaywindow_addmenubar(DisplayWindow *dw) { { "ButtonDirAxAction", "dtr-dirax", "Run _DirAx", NULL, NULL, G_CALLBACK(displaywindow_dirax) }, { "ButtonTiltAxisAction", "dtr-tiltaxis", "Tilt Axis", NULL, NULL, G_CALLBACK(displaywindow_setaxis) }, - { "ButtonRefineStepAction", "dtr-refine", "Refine Image", NULL, NULL, G_CALLBACK(displaywindow_refinestep) }, - { "ButtonRefineSeqAction", "dtr-refineall", "Refine Stack", NULL, NULL, G_CALLBACK(displaywindow_refinestack) }, + { "ButtonRefineStepAction", "dtr-refine", "Refine", NULL, NULL, G_CALLBACK(displaywindow_refinestep) }, + { "ButtonRefineSeqAction", "dtr-refineall", "Sequence", NULL, NULL, G_CALLBACK(displaywindow_refinestack) }, { "ButtonFirstImageAction", GTK_STOCK_GOTO_FIRST, "First Image", NULL, NULL, G_CALLBACK(displaywindow_image_first) }, { "ButtonPrevImageAction", GTK_STOCK_GO_BACK, "Previous Image", NULL, NULL, G_CALLBACK(displaywindow_image_prev) }, { "ButtonNextImageAction", GTK_STOCK_GO_FORWARD, "Next Image", NULL, NULL, G_CALLBACK(displaywindow_image_next) }, diff --git a/src/image.c b/src/image.c index ebeff15..d8bba0c 100644 --- a/src/image.c +++ b/src/image.c @@ -3,7 +3,7 @@ * * Handle images and image features * - * (c) 2007 Thomas White <taw27@cam.ac.uk> + * (c) 2007-2008 Thomas White <taw27@cam.ac.uk> * * dtr - Diffraction Tomography Reconstruction * diff --git a/src/image.h b/src/image.h index 6b2f4e7..62a4d2c 100644 --- a/src/image.h +++ b/src/image.h @@ -3,7 +3,7 @@ * * Handle images and image features * - * (c) 2007 Thomas White <taw27@cam.ac.uk> + * (c) 2007-2008 Thomas White <taw27@cam.ac.uk> * * dtr - Diffraction Tomography Reconstruction * diff --git a/src/iprtest.c b/src/iprtest.c deleted file mode 100644 index 369fd51..0000000 --- a/src/iprtest.c +++ /dev/null @@ -1,180 +0,0 @@ -/* - * iprtest.c - * - * Unit test for IPR - * - * (c) 2007-2008 Thomas White <taw27@cam.ac.uk> - * - * dtr - Diffraction Tomography Reconstruction - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <gtk/gtk.h> -#include <stdlib.h> -#include <math.h> - -#include "basis.h" -#include "reflections.h" -#include "image.h" -#include "reproject.h" -#include "refine.h" -#include "utils.h" -#include "mapping.h" - -typedef struct { - - ImageRecord *image; - - Basis *cell; - ReflectionList *reflections; - ReflectionList *gen; - -} IPRTestContext; - -static gint main_window_closed(GtkWidget *window, IPRTestContext *ctx) { - gtk_exit(0); - return 0; -} - -static gboolean main_window_expose(GtkWidget *drawing, GdkEventExpose *event, IPRTestContext *ctx) { - - cairo_t *cr; - ImageFeature *f; - size_t i; - - cr = gdk_cairo_create(drawing->window); - - cairo_new_path(cr); - cairo_rectangle(cr, 0.0, 0.0, 640.0, 640.0); - cairo_set_source_rgb(cr, 0.0, 0.0, 0.0); - cairo_fill(cr); - - f = ctx->image->features->features; - cairo_set_source_rgb(cr, 1.0, 1.0, 1.0); - for ( i=0; i<ctx->image->features->n_features; i++ ) { - cairo_new_path(cr); - cairo_arc(cr, f[i].x, f[i].y, 4.0, 0, 2*M_PI); - cairo_fill(cr); - } - - f = ctx->image->rflist->features; - cairo_set_line_width(cr, 1.0); - for ( i=0; i<ctx->image->rflist->n_features; i++ ) { - - cairo_new_path(cr); - cairo_arc(cr, f[i].x, f[i].y, 8.0, 0, 2*M_PI); - - if ( f[i].partner != NULL ) { - - cairo_set_source_rgb(cr, 0.0, 0.8, 0.0); - cairo_stroke(cr); - - cairo_new_path(cr); - cairo_move_to(cr, f[i].x, f[i].y); - cairo_line_to(cr, f[i].partner->x, f[i].partner->y); - cairo_stroke(cr); - - } else { - - cairo_set_source_rgb(cr, 1.0, 1.0, 1.0); - cairo_stroke(cr); - - } - - } - - cairo_destroy(cr); - - return FALSE; - -} - -static gboolean main_button_press(GtkWidget *drawing, GdkEventButton *event, IPRTestContext *ctx) { - - if ( (event->type == GDK_BUTTON_PRESS) && (event->button == 1) ) { - refine_fit_image(ctx->cell, ctx->image); - image_feature_list_free(ctx->image->rflist); - reflection_list_from_new_cell(ctx->gen, ctx->cell); - ctx->image->rflist = reproject_get_reflections(ctx->image, ctx->gen); - gdk_window_invalidate_rect(drawing->window, &drawing->allocation, FALSE); - } - - return FALSE; -} - -int main(int argc, char *argv[]) { - - GtkWidget *window; - GtkWidget *drawing; - IPRTestContext *ctx; - - gtk_init(&argc, &argv); - - ctx = malloc(sizeof(IPRTestContext)); - - ctx->image = malloc(sizeof(ImageRecord)); - ctx->image->image = NULL; - ctx->image->tilt = deg2rad(0.0); - ctx->image->omega = deg2rad(0.0); - ctx->image->slop = deg2rad(0.0); - ctx->image->fmode = FORMULATION_PIXELSIZE; - ctx->image->pixel_size = 0.6e8; - ctx->image->camera_len = 0.0; - ctx->image->lambda = lambda(1000000); - ctx->image->resolution = 1.0; - ctx->image->width = 640; - ctx->image->height = 640; - ctx->image->x_centre = 320; - ctx->image->y_centre = 320; - ctx->image->features = NULL; - ctx->image->rflist = NULL; - - ctx->cell = malloc(sizeof(Basis)); - /* The "true" cell */ - ctx->cell->a.x = 5.0e9; ctx->cell->a.y = 0.0e9; ctx->cell->a.z = 0.0e9; - ctx->cell->b.x = 0.0e9; ctx->cell->b.y = 5.0e9; ctx->cell->b.z = 0.0e9; - ctx->cell->c.x = 0.0e9; ctx->cell->c.y = 0.0e9; ctx->cell->c.z = 5.0e9; - ctx->reflections = reflection_list_from_cell(ctx->cell); - ctx->image->features = reproject_get_reflections(ctx->image, ctx->reflections); - /* The "model" cell */ - ctx->cell->a.x = 5.2e9; ctx->cell->a.y = 0.1e9; ctx->cell->a.z = 0.0e9; - ctx->cell->b.x = 0.2e9; ctx->cell->b.y = 4.8e9; ctx->cell->b.z = 0.0e9; - ctx->cell->c.x = 0.0e9; ctx->cell->c.y = 0.0e9; ctx->cell->c.z = 5.0e9; - ctx->gen = reflection_list_from_cell(ctx->cell); - ctx->image->rflist = reproject_get_reflections(ctx->image, ctx->gen); - - window = gtk_window_new(GTK_WINDOW_TOPLEVEL); - gtk_window_set_title(GTK_WINDOW(window), "IPR Test"); - gtk_widget_set_size_request(GTK_WIDGET(window), 640, 640); - g_signal_connect(G_OBJECT(window), "destroy", G_CALLBACK(main_window_closed), ctx); - gtk_window_set_resizable(GTK_WINDOW(window), FALSE); - - drawing = gtk_drawing_area_new(); - gtk_container_add(GTK_CONTAINER(window), drawing); - g_signal_connect(G_OBJECT(drawing), "expose_event", G_CALLBACK(main_window_expose), ctx); - gtk_widget_add_events(drawing, GDK_BUTTON_PRESS_MASK); - g_signal_connect(GTK_OBJECT(drawing), "button_press_event", G_CALLBACK(main_button_press), ctx); - - gtk_widget_show_all(window); - - gtk_main(); - - return 0; - -} - -/* Dummy function stubs */ -#include "displaywindow.h" -#include "gtk-valuegraph.h" -void displaywindow_update_imagestack(DisplayWindow *dw) { }; -void displaywindow_enable_cell_functions(DisplayWindow *dw, gboolean g) { }; -void displaywindow_update(DisplayWindow *dw) { }; -void displaywindow_error(const char *msg, DisplayWindow *dw) { }; -guint gtk_value_graph_get_type() { return 0; }; -GtkWidget *gtk_value_graph_new() { return NULL; }; -void gtk_value_graph_set_data(GtkValueGraph *vg, double *data, unsigned int n) { }; - 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 <stdio.h> #include <assert.h> #include <string.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_linalg.h> +#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; row<M->size1; 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; i<ctx->images->n_images; i++ ) { + int j; - for ( j=0; j<M->size2; 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; col<M->size2; col++ ) { - - int zero2 = 1; - - if ( row_altered ) break; /* This row is now fixed, so move on to the next */ - - for ( i=0; i<M->size1; 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; j<ctx->images->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; i<ctx->images->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; j<image->rflist->n_features; j++ ) { + for ( j=0; j<ctx->images->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; i<ctx->images->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; j<nf; j++ ) { - n = 0; - for ( j=0; j<image->rflist->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; i<ctx->images->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; i<ctx->images->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); diff --git a/src/refine.h b/src/refine.h index b15e56e..d3ca6d9 100644 --- a/src/refine.h +++ b/src/refine.h @@ -3,7 +3,7 @@ * * Refine the reconstruction * - * (c) 2007 Thomas White <taw27@cam.ac.uk> + * (c) 2007-2008 Thomas White <taw27@cam.ac.uk> * * dtr - Diffraction Tomography Reconstruction * @@ -17,12 +17,9 @@ #endif #include "control.h" -#include "image.h" -#include "reflections.h" -extern ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image); -extern void refine_do_stack(ControlContext *ctx); -extern void refine_do_image(ControlContext *ctx); +extern void refine_do_sequence(ControlContext *ctx); +extern void refine_do_cell(ControlContext *ctx); #endif /* REFINE_H */ |