From 5a5e330a0477f9db2f5af85637e9190be0e47fcc Mon Sep 17 00:00:00 2001 From: taw27 Date: Thu, 10 Apr 2008 18:01:34 +0000 Subject: Refinement test framework git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@271 bf6ca9ba-c028-0410-8290-897cf20841d1 --- Makefile.am | 2 ++ src/Makefile.am | 5 +++- src/refine.c | 38 +++++++++++++----------- src/refinetest.c | 88 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/reproject.c | 2 +- src/reproject.h | 2 +- 6 files changed, 117 insertions(+), 20 deletions(-) create mode 100644 src/refinetest.c diff --git a/Makefile.am b/Makefile.am index 9ac3ebb..3f6d2bf 100644 --- a/Makefile.am +++ b/Makefile.am @@ -3,3 +3,5 @@ EXTRA_DIST = configure src/displaywindow.h src/trackball.h src/reflections.h src src/itrans-zaefferer.h src/itrans-stat.h src/mapping.h src/reproject.h src/prealign.h \ src/dirax.h src/image.h src/refine.h src/gtk-valuegraph.h src/intensities.h src/glbits.h src/gtk-valuegraph.h SUBDIRS = src data +TESTS = src/refinetest + diff --git a/src/Makefile.am b/src/Makefile.am index e9aae1c..ea94888 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,10 +1,13 @@ -bin_PROGRAMS = dtr +bin_PROGRAMS = dtr refinetest 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 +refinetest_SOURCES = refinetest.c reflections.c basis.c utils.c reproject.c image.c refine.c mapping.c control.c +refinetest_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/refine.c b/src/refine.c index b56cc05..7bb2585 100644 --- a/src/refine.c +++ b/src/refine.c @@ -93,11 +93,11 @@ void refine_do_sequence(ControlContext *ctx) { gtk_widget_show_all(window_fit); /* Perform final refinement */ - refine_do_cell(ctx); printf("Best omega offset = %f deg (%f nm^-1)\n", rad2deg(omega_offs_best), fit_best/1e9); for ( j=0; jimages->n_images; j++ ) { ctx->images->images[j].omega += omega_offs_best; } + refine_do_cell(ctx); reproject_lattice_changed(ctx); mapping_adjust_axis(ctx, omega_offs_best); @@ -165,13 +165,11 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac) } -static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) { +static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug) { int v_worst, v_best, v_second_worst, i; double fom_worst, fom_new, fom_best, fom_second_worst; - -// printf("---\n"); - + /* Find the least favourable vertex of the simplex */ v_worst = 0; fom_worst = 0.0; @@ -185,7 +183,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) { fom = refine_mean_dev(d, nf, s, i); - //printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9); + if ( debug ) printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9); if ( fom > fom_worst ) { v_second_worst = v_worst; fom_second_worst = fom_worst; @@ -198,27 +196,28 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) { } } + if ( debug ) printf("The worst vertex is number %i\n", v_worst); /* Reflect this vertex across the opposite face */ refine_simplex_transform(s, v_worst, -1.0); /* Is the worst vertex any better? */ fom_new = refine_mean_dev(d, nf, s, v_worst); -// printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/1e9); + if ( debug ) printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/1e9); if ( fom_new > fom_worst ) { double fom_new_new; /* It's worse than before. Contract in 1D and see if that helps. */ -// printf("Worse. Trying a 1D contraction...\n"); + if ( debug ) printf("Worse. Trying a 1D contraction...\n"); refine_simplex_transform(s, v_worst, 0.5); fom_new_new = refine_mean_dev(d, nf, s, v_worst); -// printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/1e9); + if ( debug ) printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/1e9); if ( fom_new_new > fom_second_worst ) { int i; -// printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best); + if ( debug ) printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best); for ( i=0; i<10; i++ ) { if ( i != v_best ) refine_simplex_transform(s, i, 0.5); } @@ -227,22 +226,24 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) { } else { + #if 0 /* It's better. Try to expand in this direction */ double fom_new_new; SimplexVertex save; -// printf("This is better. Trying to expand...\n"); + if ( debug ) printf("This is better. Trying to expand...\n"); save = s[v_worst]; refine_simplex_transform(s, v_worst, 2); /* Better? */ fom_new_new = refine_mean_dev(d, nf, s, v_worst); -// printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/1e9); + if ( debug ) printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/1e9); if ( fom_new_new > fom_new ) { /* "Got too confident" */ s[v_worst] = save; -// printf("Got too confident - reverting\n"); + if ( debug ) printf("Got too confident - reverting\n"); } /* else yay. */ + #endif } @@ -257,6 +258,7 @@ double refine_do_cell(ControlContext *ctx) { double delta; int i, nf, f, it, maxiter; const double tol = 0.001e9; /* Stopping condition */ + int debug = 1; if ( !ctx->cell_lattice ) { displaywindow_error("No reciprocal unit cell has been found.", ctx->dw); @@ -282,7 +284,7 @@ double refine_do_cell(ControlContext *ctx) { } } - printf("RF: There are %i partnered features in total\n", nf); + if ( debug ) printf("RF: There are %i partnered features in total\n", nf); /* Initialise the 'deviation table' */ d = malloc(nf*sizeof(Deviation)); @@ -355,13 +357,15 @@ double refine_do_cell(ControlContext *ctx) { double conv; -// printf("Simplex method iteration %i\n", it); - conv = refine_iteration(s, d, nf); + if ( debug ) printf("------------------- Simplex method iteration %i -------------------\n", it); + conv = refine_iteration(s, d, nf, debug); if ( conv < tol ) { - printf("Converged after %i iterations (%f nm^-1)\n", it, conv/1e9); + if ( debug ) printf("RF: Converged after %i iterations (%f nm^-1)\n", it, conv/1e9); break; } + } + if ( it == maxiter ) printf("RF: Did not converge.\n"); /* Apply the final values to the cell */ ctx->cell->a.x += s[0].dax; ctx->cell->b.x += s[0].dbx; ctx->cell->c.x += s[0].dcx; diff --git a/src/refinetest.c b/src/refinetest.c new file mode 100644 index 0000000..f24862c --- /dev/null +++ b/src/refinetest.c @@ -0,0 +1,88 @@ +/* + * refinectx.c + * + * Unit test for refinement procedure + * + * (c) 2007-2008 Thomas White + * + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "basis.h" +#include "reflections.h" +#include "image.h" +#include "reproject.h" +#include "refine.h" +#include "utils.h" +#include "mapping.h" +#include "control.h" +#include "displaywindow.h" + +int main(int argc, char *argv[]) { + + ControlContext *ctx; + ReflectionList *reflections_real; + Basis *cell_real; + + ctx = control_ctx_new(); + ctx->omega = deg2rad(0.0); + ctx->lambda = lambda(300.0e3); /* 300 keV */ + ctx->fmode = FORMULATION_PIXELSIZE; + ctx->x_centre = 256; + ctx->y_centre = 256; + ctx->pixel_size = 5e7; + image_add(ctx->images, NULL, 512, 512, deg2rad(0.0), ctx); + ctx->images->images[0].features = image_feature_list_new(); + + /* Fudge to avoid horrifying pointer-related death */ + ctx->dw = malloc(sizeof(DisplayWindow)); + ctx->dw->cur_image = 0; + + /* The "true" cell */ + cell_real = malloc(sizeof(Basis)); + cell_real->a.x = 5.0e9; cell_real->a.y = 0.0e9; cell_real->a.z = 0.0e9; + cell_real->b.x = 0.0e9; cell_real->b.y = 5.0e9; cell_real->b.z = 0.0e9; + cell_real->c.x = 0.0e9; cell_real->c.y = 0.0e9; cell_real->c.z = 5.0e9; + /* The "real" reflections */ + reflections_real = reflection_list_from_cell(cell_real); + ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real); + printf("RT: %i test features generated.\n", ctx->images->images[0].features->n_features); + + /* The "model" cell to be refined */ + ctx->cell = malloc(sizeof(Basis)); + 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->cell_lattice = reflection_list_from_cell(ctx->cell); + ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice); + reproject_partner_features(ctx->images->images[0].rflist, &ctx->images->images[0]); + + refine_do_cell(ctx); + image_feature_list_free(ctx->images->images[0].rflist); + reflection_list_from_new_cell(ctx->cell_lattice, ctx->cell); + ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice); + + free(ctx); + + return 0; + +} + +/* Dummy function stubs */ +#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/reproject.c b/src/reproject.c index 4036863..e1ce13e 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -3,7 +3,7 @@ * * Synthesize diffraction patterns * - * (c) 2007 Thomas White + * (c) 2007-2008 Thomas White * * dtr - Diffraction Tomography Reconstruction * diff --git a/src/reproject.h b/src/reproject.h index 26d28bb..93f010a 100644 --- a/src/reproject.h +++ b/src/reproject.h @@ -3,7 +3,7 @@ * * Synthesize diffraction patterns * - * (c) 2007 Thomas White + * (c) 2007-2008 Thomas White * * dtr - Diffraction Tomography Reconstruction * -- cgit v1.2.3