diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-11-13 18:08:20 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-11-13 18:08:20 +0000 |
commit | 10e096c39d36e42d84a1a58a29b506aa77ee1e86 (patch) | |
tree | c2d484c27f35da5f57dd630756ce8814aa894e2a | |
parent | 54b76ad76148c5cb6093ecb16ca30e92d3c7af3a (diff) |
Refinement sequencing
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@195 bf6ca9ba-c028-0410-8290-897cf20841d1
-rw-r--r-- | src/gtk-valuegraph.c | 15 | ||||
-rw-r--r-- | src/refine.c | 157 | ||||
-rw-r--r-- | src/utils.h | 2 |
3 files changed, 157 insertions, 17 deletions
diff --git a/src/gtk-valuegraph.c b/src/gtk-valuegraph.c index b5d2bf1..555cec9 100644 --- a/src/gtk-valuegraph.c +++ b/src/gtk-valuegraph.c @@ -11,12 +11,14 @@ #include <gtk/gtk.h> #include <math.h> +#include <stdlib.h> #include "gtk-valuegraph.h" static GtkObjectClass *parent_class = NULL; static void gtk_value_graph_destroy(GtkObject *gtk_value_graph) { + free(GTK_VALUE_GRAPH(gtk_value_graph)->data); parent_class->destroy(gtk_value_graph); } @@ -161,8 +163,17 @@ static gint gtk_value_graph_draw(GtkWidget *graph, GdkEventExpose *event, gpoint cairo_move_to(cr, bw_left, height-bw_bottom-1-scale*(vg->data[0]-vg->ymin)); for ( i=1; i<vg->n; i++ ) { - cairo_line_to(cr, bw_left+((double)i/vg->xmax)*(width-bw_left-bw_right), height-bw_bottom-1-scale*(vg->data[i]-vg->ymin)); - + int p; + + p = i-1; + if ( p < 1 ) p = 1; + + if ( !isnan(vg->data[i]) && !isnan(vg->data[p]) ) { + cairo_line_to(cr, bw_left+((double)i/vg->xmax)*(width-bw_left-bw_right), height-bw_bottom-1-scale*(vg->data[i]-vg->ymin)); + } else { + cairo_move_to(cr, bw_left+((double)i/vg->xmax)*(width-bw_left-bw_right), height-bw_bottom-1-scale*(vg->data[i]-vg->ymin)); + } + } cairo_set_line_width(cr, 0.5); cairo_set_source_rgb(cr, 0.0, 0.0, 1.0); diff --git a/src/refine.c b/src/refine.c index b2dd5d0..85573b6 100644 --- a/src/refine.c +++ b/src/refine.c @@ -18,6 +18,7 @@ #include <stdlib.h> #include <stdio.h> #include <assert.h> +#include <string.h> #include "displaywindow.h" #include "gtk-valuegraph.h" @@ -28,6 +29,7 @@ #include "control.h" #include "mapping.h" #include "imagedisplay.h" +#include "utils.h" /* Return the root sum squared deviation distance for all the "reprojectable" features in an image */ static double refine_image_deviation(ImageRecord *image, ReflectionList *cell_lattice) { @@ -64,6 +66,7 @@ typedef enum { INDEX_C = 1<<2 } RefinementIndex; +#if 0 static const char *refine_decode(RefinementIndex i) { switch ( i ) { @@ -81,6 +84,7 @@ static const char *refine_decode(RefinementIndex i) { } } +#endif /* Use the IPR algorithm to make "cell" fit the given image */ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *cell_lattice) { @@ -123,7 +127,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio /* Determine the difference vector */ dix = rflist->features[i].partner->x - rflist->features[i].x; diy = rflist->features[i].partner->y - rflist->features[i].y; - printf("Feature %3i: %3i %3i %3i dev = %f %f px\n", i, h, k, l, dix, diy); + // 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 = rflist->features[i].partner->x; @@ -133,7 +137,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio mapping_map_to_space(rflist->features[i].partner, &dx, &dy, &dz, &twotheta); rflist->features[i].partner->x = old_x; rflist->features[i].partner->y = old_y; - printf("dev=%8e %8e %8e (%5f mrad)\n", dx, dy, dz, twotheta*1e3); + // 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; @@ -152,7 +156,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio 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("dev(hkl) = %f %f %f\n", dh, dk, dl); + // 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; @@ -161,17 +165,17 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio shared = index & delta; - printf(" index: %s\n delta: %s\nshared: %s\n", refine_decode(index), refine_decode(delta), refine_decode(shared)); + // 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("Pure shear.\n"); + // printf("RF: Pure shear.\n"); } if ( shared & INDEX_A ) { double w = (double)abs(h) / (abs(h)+abs(k)+abs(l)); - printf("w(a) = %f\n", w); + // 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; @@ -179,7 +183,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio } if ( shared & INDEX_B ) { double w = (double)abs(k) / (abs(h)+abs(k)+abs(l)); - printf("w(b) = %f\n", w); + // 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; @@ -187,16 +191,16 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio } if ( shared & INDEX_C ) { double w = (double)abs(l) / (abs(h)+abs(k)+abs(l)); - printf("w(c) = %f\n", w); + // 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("Distortion along x: %+8e = %+8e\n", h*cd.a.x + k*cd.b.x + l*cd.c.x, dx); - printf("Distortion along y: %+8e = %+8e\n", h*cd.a.y + k*cd.b.y + l*cd.c.y, dy); - printf("Distortion along z: %+8e = %+8e\n", h*cd.a.z + k*cd.b.z + l*cd.c.z, dz); + // 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); done = TRUE; @@ -207,7 +211,7 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio return NULL; } - printf("n_a=%i, n_b=%i, n_c=%i\n", n_a, n_b, n_c); + //printf("RF: n_a=%i, n_b=%i, n_c=%i\n", n_a, n_b, n_c); if ( n_a ) { cd.a.x /= (double)n_a; cd.a.y /= (double)n_a; cd.a.z /= (double)n_a; } @@ -218,9 +222,9 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio cd.c.x /= (double)n_c; cd.c.y /= (double)n_c; cd.c.z /= (double)n_c; } - printf("Total distortion(a) = %+8e %+8e %+8e\n", cd.a.x, cd.a.y, cd.a.z); - printf("Total distortion(b) = %+8e %+8e %+8e\n", cd.b.x, cd.b.y, cd.b.z); - printf("Total distortion(c) = %+8e %+8e %+8e\n", cd.c.x, cd.c.y, cd.c.z); +// printf("RF: Total distortion(a) = %+8e %+8e %+8e\n", cd.a.x, cd.a.y, cd.a.z); +// printf("RF: Total distortion(b) = %+8e %+8e %+8e\n", cd.b.x, cd.b.y, cd.b.z); +// printf("RF: Total distortion(c) = %+8e %+8e %+8e\n", cd.c.x, cd.c.y, cd.c.z); cell->a.x += cd.a.x; cell->a.y += cd.a.y; cell->a.z += cd.a.z; cell->b.x += cd.b.x; cell->b.y += cd.b.y; cell->b.z += cd.b.z; @@ -291,14 +295,137 @@ static gint refine_step(GtkWidget *step_button, ControlContext *ctx) { } +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; + + for ( i=0; i<ctx->images->n_images; i++ ) { + + /* Ensure lattice is up to date */ + reproject_lattice_changed(ctx); + + if ( is_odd(i) ) { + + /* Odd-numbered images: measure */ + ImageFeatureList *rflist; + ImageFeatureList *flist; + ImageRecord *image; + int j, n; + double image_dev_mean = 0; + + image = &ctx->images->images[i]; + rflist = reproject_get_reflections(image, ctx->cell_lattice); + flist = image->features; + reproject_partner_features(rflist, image); + + n = 0; + for ( j=0; j<rflist->n_features; j++ ) { + + double dix, diy; + + /* Skip if no partner */ + if ( !rflist->features[j].partner ) continue; + + /* Determine the difference vector */ + dix = rflist->features[j].partner->x - rflist->features[j].x; + diy = rflist->features[j].partner->y - 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++; + + image_feature_list_free(rflist); + + } else { + + /* Even-numbered images: fit */ + refine_fit_image(ctx->cell, &ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice); + displaywindow_update(ctx->dw); + + } + + } + + series_dev_mean /= series_dev_n; + *fit = series_dev_mean; + *warp = (series_dev_max - series_dev_min)/series_dev_min; + + return 0; + +} + static gint refine_sequence(GtkWidget *step_button, 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 ) { displaywindow_error("No reciprocal unit cell has been found.", ctx->dw); return 0; } + 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"); + 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; i<ctx->images->n_images; i++ ) { + ctx->images->images[i].omega -= omega_offs; + } + memcpy(ctx->cell, &cell_copy, sizeof(Basis)); + + } + 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; diff --git a/src/utils.h b/src/utils.h index 55e98f2..ab551e9 100644 --- a/src/utils.h +++ b/src/utils.h @@ -33,5 +33,7 @@ extern void chomp(char *s); #define rad2deg(a) ((a)*180/M_PI) #define deg2rad(a) ((a)*M_PI/180) +#define is_odd(a) ((a)%2==1) + #endif /* UTILS_H */ |