diff options
-rw-r--r-- | src/displaywindow.c | 13 | ||||
-rw-r--r-- | src/displaywindow.h | 1 | ||||
-rw-r--r-- | src/intensities.c | 2 | ||||
-rw-r--r-- | src/refine.c | 634 | ||||
-rw-r--r-- | src/refinetest2d.c | 8 | ||||
-rw-r--r-- | src/refinetest3d1.c | 126 | ||||
-rw-r--r-- | src/refinetest3d2.c | 130 | ||||
-rw-r--r-- | src/refinetest3d3.c | 108 | ||||
-rw-r--r-- | src/reflections.h | 3 | ||||
-rw-r--r-- | src/reproject.h | 3 |
10 files changed, 563 insertions, 465 deletions
diff --git a/src/displaywindow.c b/src/displaywindow.c index 1ffd95e..d0c8e3e 100644 --- a/src/displaywindow.c +++ b/src/displaywindow.c @@ -372,7 +372,15 @@ static gint displaywindow_image_last(GtkWidget *widget, DisplayWindow *dw) } static gint displaywindow_refinestep(GtkWidget *widget, DisplayWindow *dw) { + + int old_n; + + old_n = dw->ctx->images->n_images; + dw->ctx->images->n_images = 1; + printf("WARNING: I'm only refining based on the first image\n"); refine_do_cell(dw->ctx); + dw->ctx->images->n_images = old_n; + return 0; } @@ -573,7 +581,8 @@ void displaywindow_update_imagestack(DisplayWindow *dw) { /* Perform relrod projection if necessary */ if ( !image->rflist ) { image->rflist = reproject_get_reflections(image, - dw->ctx->cell_lattice); + dw->ctx->cell_lattice, + dw->ctx); } /* Draw the reprojected peaks */ @@ -680,7 +689,7 @@ static void displaywindow_scrolltoend(DisplayWindow *dw) dw->messages_mark, 0, TRUE, 1.0, 0.0); } -static void displaywindow_message(DisplayWindow *dw, const char *text) +void displaywindow_message(DisplayWindow *dw, const char *text) { GtkTextBuffer *buffer; GtkTextIter iter; diff --git a/src/displaywindow.h b/src/displaywindow.h index a69e807..8330f7e 100644 --- a/src/displaywindow.h +++ b/src/displaywindow.h @@ -100,6 +100,7 @@ extern void displaywindow_update(DisplayWindow *dw); extern void displaywindow_update_dirax(ControlContext *ctx, DisplayWindow *dw); extern void displaywindow_error(const char *msg, DisplayWindow *dw); extern void displaywindow_enable_cell_functions(DisplayWindow *dw, gboolean g); +extern void displaywindow_message(DisplayWindow *dw, const char *text); #endif /* DISPLAYWINDOW_H */ diff --git a/src/intensities.c b/src/intensities.c index 2a0bd88..88de5f2 100644 --- a/src/intensities.c +++ b/src/intensities.c @@ -136,7 +136,7 @@ void intensities_extract(ControlContext *ctx) image = &ctx->images->images[i]; if ( image->rflist == NULL ) image->rflist = reproject_get_reflections(image, - ctx->cell_lattice); + ctx->cell_lattice, ctx); /* For each reprojected feature */ for ( j=0; j<image->rflist->n_features; j++ ) { diff --git a/src/refine.c b/src/refine.c index 3a61170..f3c1dc3 100644 --- a/src/refine.c +++ b/src/refine.c @@ -3,13 +3,13 @@ * * Refine the reconstruction * - * (c) 2007-2008 Thomas White <taw27@cam.ac.uk> + * (c) 2007-2009 Thomas White <taw27@cam.ac.uk> * * dtr - Diffraction Tomography Reconstruction * */ -#ifdef HAVE_CONFIG_H +#if HAVE_CONFIG_H #include <config.h> #endif @@ -36,22 +36,43 @@ #define NUM_PARAMS 9 /* Refine debug */ -#define REFINE_DEBUG 1 +#define REFINE_DEBUG_MORE 0 -/* 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; +#define DELTA (0.0001e9) typedef struct { + /* Indicies of reflection */ signed int h; signed int k; signed int l; - double dx; double dy; double dz; + /* Deviation vector for this reflection in nm^-1 */ + double dx; double dy; + /* x,y,z components of "image x" */ + double xx; double xy; double xz; + /* x,y,z components of "image y" */ + double yx; double yy; double yz; + + #if REFINE_DEBUG_MORE + /* Image number this dev came from, for debugging purposes */ + int image; + double rdx; double rdy; + #endif /* REFINE_DEBUG_MORE */ + } Deviation; -void refine_do_sequence(ControlContext *ctx) { +static void refine_printf(DisplayWindow *dw, char *fmt, ...) +{ + va_list ap; + char tmp[1024]; + + va_start(ap, fmt); + vsnprintf(tmp, 1023, fmt, ap); + va_end(ap); + printf("%s", tmp); + displaywindow_message(dw, tmp); +} + +void refine_do_sequence(ControlContext *ctx) +{ double omega_offs; int idx; double *fit_vals; @@ -59,41 +80,43 @@ void refine_do_sequence(ControlContext *ctx) { GtkWidget *graph_fit; double fit_best, omega_offs_best; int j; - + fit_vals = malloc(401*sizeof(double)); idx = 0; - + fit_best = 1000.0e9; omega_offs_best = 0.0; - for ( omega_offs=-deg2rad(2.0); omega_offs<=deg2rad(2.0); omega_offs+=deg2rad(0.01) ) { + for ( omega_offs=-deg2rad(2.0); omega_offs<=deg2rad(2.0); + omega_offs+=deg2rad(0.01) ) { double fit; int i; Basis cell_copy; - + cell_copy = *ctx->cell; - + for ( i=0; i<ctx->images->n_images; i++ ) { ctx->images->images[i].omega += omega_offs; } reproject_lattice_changed(ctx); - + fit = refine_do_cell(ctx); - printf("RF: omega_offs=%f deg, fit=%f nm^-1\n", rad2deg(omega_offs), - fit/DISPFACTOR); + refine_printf(ctx->dw, " " + "omega_offs=%f deg, fit=%f nm^-1\n", + rad2deg(omega_offs), fit/DISPFACTOR); fit_vals[idx++] = fit; if ( fit < fit_best ) { fit_best = fit; omega_offs_best = omega_offs; } - + for ( i=0; i<ctx->images->n_images; i++ ) { ctx->images->images[i].omega -= omega_offs; } *ctx->cell = cell_copy; } - + 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"); @@ -101,367 +124,376 @@ void refine_do_sequence(ControlContext *ctx) { 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); - + /* Perform final refinement */ - printf("Best omega offset = %f deg (%f nm^-1)\n", rad2deg(omega_offs_best), fit_best/DISPFACTOR); + refine_printf(ctx->dw, "Best omega offset = %f deg (%f nm^-1)\n", + rad2deg(omega_offs_best), fit_best/DISPFACTOR); for ( j=0; j<ctx->images->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); - } -static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) { - +static double refine_mean_dev(Deviation *d, int nf, Basis *devcell, int disp) +{ double fom = 0.0; int f; - + for ( f=0; f<nf; f++ ) { - + double xdf, ydf, zdf; - - xdf = d[f].h*s[i].dax + d[f].k*s[i].dbx + d[f].l*s[i].dcx; - ydf = d[f].h*s[i].day + d[f].k*s[i].dby + d[f].l*s[i].dcy; - zdf = d[f].h*s[i].daz + d[f].k*s[i].dbz + d[f].l*s[i].dcz; - xdf -= d[f].dx; - ydf -= d[f].dy; - zdf -= d[f].dz; - - fom += sqrt(xdf*xdf + ydf*ydf + zdf*zdf); - - } - -// return fabs(s[i].dax-10.0) + fabs(s[i].day-20.0) + fabs(s[i].dbx-30.0) + fabs(s[i].dby-40.0); - - return fom/nf; + double dx, dy; -} + /* Calculate 3D deviation vector */ + xdf = d[f].h*devcell->a.x + d[f].k*devcell->b.x + + d[f].l*devcell->c.x; + ydf = d[f].h*devcell->a.y + d[f].k*devcell->b.y + + d[f].l*devcell->c.y; + zdf = d[f].h*devcell->a.z + d[f].k*devcell->b.z + + d[f].l*devcell->c.z; -static void refine_display_simplex(SimplexVertex sv) { + /* Project into original image */ + dx = xdf*d[f].xx + ydf*d[f].xy + zdf*d[f].xz; + dy = xdf*d[f].yx + ydf*d[f].yy + zdf*d[f].yz; - printf("%9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", - sv.dax/DISPFACTOR, sv.day/DISPFACTOR, sv.daz/DISPFACTOR, - sv.dbx/DISPFACTOR, sv.dby/DISPFACTOR, sv.dbz/DISPFACTOR, - sv.dcx/DISPFACTOR, sv.dcy/DISPFACTOR, sv.dcz/DISPFACTOR); - -} + dx -= d[f].dx; + dy -= d[f].dy; -/* Expand the simplex across from vertex v_worst by factor 'fac'. - * fac = -1 is a reflection - * fac = +n is a 1d expansion - */ -static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac) { - - SimplexVertex centre; - int i, nv; - - /* Average the coordinates of the non-worst vertices to - * get the centre of the opposite face. */ - centre.dax = 0.0; centre.day = 0.0; centre.daz = 0.0; - centre.dbx = 0.0; centre.dby = 0.0; centre.dbz = 0.0; - centre.dcx = 0.0; centre.dcy = 0.0; centre.dcz = 0.0; - nv = 0; - for ( i=0; i<=NUM_PARAMS; i++ ) { - if ( i != v_worst ) { - centre.dax += s[i].dax; centre.day += s[i].day; centre.daz += s[i].daz; - centre.dbx += s[i].dbx; centre.dby += s[i].dby; centre.dbz += s[i].dbz; - centre.dcx += s[i].dcx; centre.dcy += s[i].dcy; centre.dcz += s[i].dcz; - nv++; + #if REFINE_DEBUG_MORE + if ( disp ) { + d[f].rdx = dx; + d[f].rdy = dy; + printf("Image %2i ref %3i %3i %3i dev %8.5f %8.5f\n", + d[f].image, d[f].h, d[f].k, d[f].l, + dx/DISPFACTOR, dy/DISPFACTOR); } + #endif /* REFINE_DEBUG_MORE */ + + fom += sqrt(dx*dx + dy*dy); + } - centre.dax /= nv; centre.day /= nv; centre.daz /= nv; - centre.dbx /= nv; centre.dby /= nv; centre.dbz /= nv; - centre.dcx /= nv; centre.dcy /= nv; centre.dcz /= nv; - -// printf("Before transformation: "); -// refine_display_simplex(s[v_worst]); - -// printf(" Midpoint: "); -// refine_display_simplex(centre); - - /* Do the transformation */ - s[v_worst].dax = centre.dax + fac * (s[v_worst].dax - centre.dax); - s[v_worst].day = centre.day + fac * (s[v_worst].day - centre.day); - s[v_worst].daz = centre.daz + fac * (s[v_worst].daz - centre.daz); - s[v_worst].dbx = centre.dbx + fac * (s[v_worst].dbx - centre.dbx); - s[v_worst].dby = centre.dby + fac * (s[v_worst].dby - centre.dby); - s[v_worst].dbz = centre.dbz + fac * (s[v_worst].dbz - centre.dbz); - s[v_worst].dcx = centre.dcx + fac * (s[v_worst].dcx - centre.dcx); - s[v_worst].dcy = centre.dcy + fac * (s[v_worst].dcy - centre.dcy); - s[v_worst].dcz = centre.dcz + fac * (s[v_worst].dcz - centre.dcz); - - if ( REFINE_DEBUG ) { - printf(" After transformation: "); - refine_display_simplex(s[v_worst]); + + return fom/nf; + +} + +static void refine_cell_delta(Basis *devcell_try, int comp) { + + switch ( comp ) { + case 0 : break; + case 1 : devcell_try->a.x += DELTA; break; + case 2 : devcell_try->a.y += DELTA; break; + case 3 : devcell_try->a.z += DELTA; break; + case 4 : devcell_try->b.x += DELTA; break; + case 5 : devcell_try->b.y += DELTA; break; + case 6 : devcell_try->b.z += DELTA; break; + case 7 : devcell_try->c.x += DELTA; break; + case 8 : devcell_try->c.y += DELTA; break; + case 9 : devcell_try->c.z += DELTA; break; + case 10 : devcell_try->a.x -= DELTA; break; + case 11 : devcell_try->a.y -= DELTA; break; + case 12 : devcell_try->a.z -= DELTA; break; + case 13 : devcell_try->b.x -= DELTA; break; + case 14 : devcell_try->b.y -= DELTA; break; + case 15 : devcell_try->b.z -= DELTA; break; + case 16 : devcell_try->c.x -= DELTA; break; + case 17 : devcell_try->c.y -= DELTA; break; + case 18 : devcell_try->c.z -= DELTA; break; + default : fprintf(stderr, + "refine_cell_delta: argument out of range\n"); } - + } -/* Contract vertex v of simplex s towards vertex v_best */ -static void refine_simplex_contract(SimplexVertex *s, int v, int v_best) { - - s[v].dax = s[v_best].dax + 0.5 * (s[v].dax - s[v_best].dax); - s[v].day = s[v_best].day + 0.5 * (s[v].day - s[v_best].day); - s[v].daz = s[v_best].daz + 0.5 * (s[v].daz - s[v_best].daz); - s[v].dbx = s[v_best].dbx + 0.5 * (s[v].dbx - s[v_best].dbx); - s[v].dby = s[v_best].dby + 0.5 * (s[v].dby - s[v_best].dby); - s[v].dbz = s[v_best].dbz + 0.5 * (s[v].dbz - s[v_best].dbz); - s[v].dcx = s[v_best].dcx + 0.5 * (s[v].dcx - s[v_best].dcx); - s[v].dcy = s[v_best].dcy + 0.5 * (s[v].dcy - s[v_best].dcy); - s[v].dcz = s[v_best].dcz + 0.5 * (s[v].dcz - s[v_best].dcz); - +static void refine_show_cell(DisplayWindow *dw, Basis cell) { + + refine_printf(dw, "a: %+10.8f %+10.8f %+10.8f\n", + cell.a.x/DISPFACTOR, cell.a.y/DISPFACTOR, + cell.a.z/DISPFACTOR); + refine_printf(dw, "b: %+10.8f %+10.8f %+10.8f\n", + cell.b.x/DISPFACTOR, cell.b.y/DISPFACTOR, + cell.b.z/DISPFACTOR); + refine_printf(dw, "c: %+10.8f %+10.8f %+10.8f\n", + cell.c.x/DISPFACTOR, cell.c.y/DISPFACTOR, + cell.c.z/DISPFACTOR); + } -static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) { - - int v_worst, v_best, v_second_worst, i; - double fom_worst, fom_new, fom_best, fom_second_worst; - - /* Find the least favourable vertex of the simplex */ - v_worst = 0; - fom_worst = 0.0; - v_best = 0; - fom_best = 100e9; - v_second_worst = 0; - fom_second_worst = 0.0; - if ( REFINE_DEBUG ) printf("Vertex FoM/nm^-1\n"); - for ( i=0; i<=NUM_PARAMS; i++ ) { - - double fom; - - fom = refine_mean_dev(d, nf, s, i); - - if ( REFINE_DEBUG ) printf("%6i %8f\n", i, fom/DISPFACTOR); - if ( fom > fom_worst ) { - v_second_worst = v_worst; - fom_second_worst = fom_worst; - fom_worst = fom; - v_worst = i; - } - if ( fom < fom_best ) { - fom_best = fom; - v_best = i; - } - - } - if ( REFINE_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); - if ( REFINE_DEBUG ) printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", - fom_new/DISPFACTOR); - if ( fom_new > fom_worst ) { - - double fom_new_new; - - /* It's worse than before. Contract in 1D and see if that helps. */ - if ( REFINE_DEBUG ) printf("Worse. Trying a 1D contraction...\n"); - /* Minus puts it back on the original side of the 'good' face */ - refine_simplex_transform(s, v_worst, -0.5); - fom_new_new = refine_mean_dev(d, nf, s, v_worst); - if ( REFINE_DEBUG ) printf("Mean deviation after 1D contraction is %f nm^-1\n", - fom_new_new/DISPFACTOR); - if ( fom_new_new > fom_second_worst ) { - - int i; - - if ( REFINE_DEBUG ) printf("Not as good as the second worst vertex: contracting around the " - "best vertex (%i)\n", v_best); - for ( i=0; i<=NUM_PARAMS; i++ ) { - if ( i != v_best ) refine_simplex_contract(s, i, v_best); - } - - } - - } else if ( fom_new < fom_worst ) { - - /* It's better. Try to expand in this direction */ - double fom_new_new; - SimplexVertex save; - - if ( REFINE_DEBUG ) printf("This is better. Trying to expand...\n"); - - save = s[v_worst]; - refine_simplex_transform(s, v_worst, 2.0); /* +ve means stay on this side of the 'good' face */ - /* Better? */ - fom_new_new = refine_mean_dev(d, nf, s, v_worst); - if ( REFINE_DEBUG ) printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/DISPFACTOR); - if ( fom_new_new > fom_new ) { - /* "Got too confident" */ - s[v_worst] = save; - if ( REFINE_DEBUG ) printf("Got too confident - reverting\n"); - } else { - if ( REFINE_DEBUG ) printf("Better still. Great.\n"); - } - - } else { - - printf("No change!\n"); - - } - - /* Check convergence and return */ - fom_worst = 0.0; - fom_best = 100e9; - for ( i=0; i<=NUM_PARAMS; i++ ) { - double fom; - fom = refine_mean_dev(d, nf, s, i); - if ( fom > fom_worst ) { - fom_worst = fom; - } - if ( fom < fom_best ) { - fom_best = fom; - } - } - - printf("Simplex size: %e - %e = %e\n", fom_worst/DISPFACTOR, fom_best/DISPFACTOR, (fom_worst - fom_best)/DISPFACTOR); - - return fom_best; +static Basis refine_test_confidence(Deviation *d, int nf) { + + Basis *devcell; + double val, new_val; + Basis confidence; + + devcell = malloc(sizeof(Basis)); + devcell->a.x = 0.0e9; devcell->b.x = 0.0e9; devcell->c.x = 0.0e9; + devcell->a.y = 0.0e9; devcell->b.y = 0.0e9; devcell->c.y = 0.0e9; + devcell->a.z = 0.0e9; devcell->b.z = 0.0e9; devcell->c.z = 0.0e9; + val = refine_mean_dev(d, nf, devcell, 0); + + devcell->a.x = 1.0e9; + new_val = refine_mean_dev(d, nf, devcell, 0); + confidence.a.x = new_val-val; + devcell->a.x = 0.0; + + devcell->a.y = 1.0e9; + new_val = refine_mean_dev(d, nf, devcell, 0); + confidence.a.y = new_val-val; + devcell->a.y = 0.0; + + devcell->a.z = 1.0e9; + new_val = refine_mean_dev(d, nf, devcell, 0); + confidence.a.z = new_val-val; + devcell->a.z = 0.0; + + devcell->b.x = 1.0e9; + new_val = refine_mean_dev(d, nf, devcell, 0); + confidence.b.x = new_val-val; + devcell->b.x = 0.0; + + devcell->b.y = 1.0e9; + new_val = refine_mean_dev(d, nf, devcell, 0); + confidence.b.y = new_val-val; + devcell->b.y = 0.0; + + devcell->b.z = 1.0e9; + new_val = refine_mean_dev(d, nf, devcell, 0); + confidence.b.z = new_val-val; + devcell->b.z = 0.0; + + devcell->c.x = 1.0e9; + new_val = refine_mean_dev(d, nf, devcell, 0); + confidence.c.x = new_val-val; + devcell->c.x = 0.0; + + devcell->c.y = 1.0e9; + new_val = refine_mean_dev(d, nf, devcell, 0); + confidence.c.y = new_val-val; + devcell->c.y = 0.0; + + devcell->c.z = 1.0e9; + new_val = refine_mean_dev(d, nf, devcell, 0); + confidence.c.z = new_val-val; + devcell->c.z = 0.0; + + return confidence; } double refine_do_cell(ControlContext *ctx) { - SimplexVertex s[10]; Deviation *d; - double delta; - int i, nf, f, it, maxiter; - const double tol = 0.001e9; /* Stopping condition */ - //const double tol = 0.001; /* For testing */ - + int i, nf, f; + Basis *devcell; + Basis *devcell_try; + Basis *devcell_try_best; + double mean_dev; + double mean_dev_try; + int it; + Basis conf; + double conf_threshold; + if ( !ctx->cell_lattice ) { - displaywindow_error("No reciprocal unit cell has been found.", ctx->dw); + displaywindow_error("No reciprocal unit cell has been found.", + ctx->dw); return -1; } - + if ( ctx->images->n_images == 0 ) { - displaywindow_error("There are no images to refine against.", ctx->dw); + displaywindow_error("There are no images to refine against.", + ctx->dw); return -1; } - + /* Determine the size of the 'deviation table' */ nf = 0; for ( i=0; i<ctx->images->n_images; i++ ) { - + int j; - + if ( !ctx->images->images[i].rflist ) { - ctx->images->images[i].rflist = reproject_get_reflections(&ctx->images->images[i], - ctx->cell_lattice); + ctx->images->images[i].rflist = + reproject_get_reflections(&ctx->images->images[i], + ctx->cell_lattice, ctx); } - + for ( j=0; j<ctx->images->images[i].rflist->n_features; j++ ) { - if ( ctx->images->images[i].rflist->features[j].partner != NULL ) nf++; + if ( ctx->images->images[i].rflist->features[j].partner + != NULL ) nf++; } - - printf("%i features from image %i\n", nf, i); - + + refine_printf(ctx->dw, "%i features from image %i\n", nf, i); + } - if ( REFINE_DEBUG ) printf("RF: There are %i partnered features in total\n", nf); - + refine_printf(ctx->dw, "There are %i partnered features in total\n", + nf); + /* Initialise the 'deviation table' */ d = malloc(nf*sizeof(Deviation)); f = 0; for ( i=0; i<ctx->images->n_images; i++ ) { - + ImageRecord *image; int j; - + image = &ctx->images->images[i]; - + for ( j=0; j<ctx->images->images[i].rflist->n_features; j++ ) { - + ImageFeature *rf; double dix, diy, dx, dy; - double dlx, dly, dlz; double old_x, old_y; - + rf = &image->rflist->features[j]; if ( !rf->partner ) continue; - + 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, d[f].h, d[f].k, d[f].l, - dix, diy); - + refine_printf(ctx->dw, + "%3i: %3i: %3i %3i %3i dev %+9.5f %+9.5f px ", + i, 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; double mod = sqrt(dx*dx + dy*dy)/DISPFACTOR; - printf("=> %+10.5f %+10.5f %+10.5f nm^-1 (length %10.5f nm^1)\n", dlx/DISPFACTOR, dly/DISPFACTOR, dlz/DISPFACTOR, mod); - - d[f].dx = dlx; - d[f].dy = dly; - d[f].dz = dlz; - + refine_printf(ctx->dw, + "= %+10.5f %+10.5f nm^-1 (length %7.5f nm^1)\n", + dx/DISPFACTOR, dy/DISPFACTOR, mod); + + d[f].dx = dx; + d[f].dy = dy; + + /* Store plane normal for use in projection later */ + mapping_rotate(1.0, 0.0, 0.0, + &d[f].xx, &d[f].xy, &d[f].xz, + image->omega, image->tilt); + mapping_rotate(0.0, 1.0, 0.0, + &d[f].yx, &d[f].yy, &d[f].yz, + image->omega, image->tilt); + + #if REFINE_DEBUG_MORE + d[f].image = i; + #endif /* REFINE_DEBUG_MORE */ + f++; - + } - + } assert( f == nf ); - - /* Initialise the simplex */ - delta = 0.01e9; - 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].day = delta; - memcpy(&s[3], &s[0], sizeof(SimplexVertex)); s[3].dbx = delta; - memcpy(&s[4], &s[0], sizeof(SimplexVertex)); s[4].dby = delta; /* 2d vertices first */ - memcpy(&s[5], &s[0], sizeof(SimplexVertex)); s[5].daz = delta; - memcpy(&s[6], &s[0], sizeof(SimplexVertex)); s[6].dbz = delta; - memcpy(&s[7], &s[0], sizeof(SimplexVertex)); s[7].dcx = delta; - memcpy(&s[8], &s[0], sizeof(SimplexVertex)); s[8].dcy = delta; - memcpy(&s[9], &s[0], sizeof(SimplexVertex)); s[9].dcz = delta; - - /* Iterate */ - maxiter = 10000; - for ( it=0; it<maxiter; it++ ) { - - double conv; - - //for ( i=0; i<5; i++ ) { - // refine_display_simplex(s[i]); - //} - - if ( REFINE_DEBUG ) printf("------------------- Simplex method iteration %i -------------------\n", it); - conv = refine_iteration(s, d, nf); - if ( conv < tol ) { - if ( REFINE_DEBUG ) printf("RF: Converged after %i iterations (%f nm^-1)\n", it, - conv/DISPFACTOR); + + /* Initial situation */ + devcell = malloc(sizeof(Basis)); + devcell->a.x = 0.0e9; devcell->b.x = 0.0e9; devcell->c.x = 0.0e9; + devcell->a.y = 0.0e9; devcell->b.y = 0.0e9; devcell->c.y = 0.0e9; + devcell->a.z = 0.0e9; devcell->b.z = 0.0e9; devcell->c.z = 0.0e9; + mean_dev = refine_mean_dev(d, nf, devcell, 0); + refine_printf(ctx->dw, "Initial mean deviation: %13.8f nm^1\n", + mean_dev/DISPFACTOR); + + /* Test confidence */ + conf = refine_test_confidence(d, nf); + + /* Determine direction of steepest gradient */ + devcell_try = malloc(sizeof(Basis)); + devcell_try_best = malloc(sizeof(Basis)); + for ( it=1; it<=10000; it++ ) { + + int found = 0; + int comp1, comp2, comp3; + + #if REFINE_DEBUG_MORE + refine_mean_dev(d, nf, devcell, 1); + printf(ctx->dw, "Iteration %i starts at dev %14.12f\n", + it, mean_dev/DISPFACTOR); + printf(ctx->dw, "Current dev cell:\n"); + printf(ctx->dw, "a = %8.5f %8.5f %8.5f\n", + devcell->a.x/DISPFACTOR, + devcell->a.y/DISPFACTOR, + devcell->a.z/DISPFACTOR); + printf(ctx->dw, "b = %8.5f %8.5f %8.5f\n", + devcell->b.x/DISPFACTOR, + devcell->b.y/DISPFACTOR, + devcell->b.z/DISPFACTOR); + printf(ctx->dw, "c = %8.5f %8.5f %8.5f\n", + devcell->c.x/DISPFACTOR, + devcell->c.y/DISPFACTOR, + devcell->c.z/DISPFACTOR); + #endif /* REFINE_DEBUG_MORE */ + + for ( comp1=1; comp1<19; comp1++ ) { + for ( comp2=0; comp2<19; comp2++ ) { + for ( comp3=0; comp3<19; comp3++ ) { + + memcpy(devcell_try, devcell, sizeof(Basis)); + refine_cell_delta(devcell_try, comp1); + refine_cell_delta(devcell_try, comp2); + refine_cell_delta(devcell_try, comp3); + mean_dev_try = refine_mean_dev(d, nf, devcell_try, 0); + + /* Improvement greater than the tolerance? */ + if ( mean_dev_try < mean_dev-0.001 ) { + mean_dev = mean_dev_try; + memcpy(devcell_try_best, devcell_try, + sizeof(Basis)); + found = 1; + } + + } + } + } +// printf("mean_dev = %f\n", mean_dev); + if ( found ) { + memcpy(devcell, devcell_try_best, sizeof(Basis)); + } else { + printf("No further change after %i iterations\n", it); break; } - + if ( !(it % 1000) ) { + refine_printf(ctx->dw, + "After %5i iterations: mean dev = %13.8f nm^1\n", + it, mean_dev/DISPFACTOR); + } + } - if ( it == maxiter ) printf("RF: Did not converge.\n"); - + free(devcell_try); + free(devcell_try_best); + refine_printf(ctx->dw, + "Final mean dev (%5i iterations) = %13.8f nm^1\n", + it, mean_dev/DISPFACTOR); + refine_printf(ctx->dw, "Final cell deviation:\n"); + refine_show_cell(ctx->dw, *devcell); + refine_printf(ctx->dw, "Confidence:\n"); + refine_show_cell(ctx->dw, conf); + /* 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; - ctx->cell->a.y += s[0].day; ctx->cell->b.y += s[0].dby; ctx->cell->c.y += s[0].dcy; - ctx->cell->a.z += s[0].daz; ctx->cell->b.z += s[0].dbz; ctx->cell->c.z += s[0].dcz; - + conf_threshold = 2.5; + if ( conf.a.x > conf_threshold ) ctx->cell->a.x += devcell->a.x; + if ( conf.a.y > conf_threshold ) ctx->cell->a.y += devcell->a.y; + if ( conf.a.z > conf_threshold ) ctx->cell->a.z += devcell->a.z; + if ( conf.b.x > conf_threshold ) ctx->cell->b.x += devcell->b.x; + if ( conf.b.y > conf_threshold ) ctx->cell->b.y += devcell->b.y; + if ( conf.b.z > conf_threshold ) ctx->cell->b.z += devcell->b.z; + if ( conf.c.x > conf_threshold ) ctx->cell->c.x += devcell->c.x; + if ( conf.c.y > conf_threshold ) ctx->cell->c.y += devcell->c.y; + if ( conf.c.z > conf_threshold ) ctx->cell->c.z += devcell->c.z; + ctx->images->images[ctx->dw->cur_image].rflist = NULL; reproject_lattice_changed(ctx); displaywindow_update(ctx->dw); - - return refine_mean_dev(d, nf, s, 0); - -} + return mean_dev; + +} diff --git a/src/refinetest2d.c b/src/refinetest2d.c index 1c30348..2b09ef4 100644 --- a/src/refinetest2d.c +++ b/src/refinetest2d.c @@ -107,6 +107,7 @@ int main(int argc, char *argv[]) { ctx->x_centre = 256; ctx->y_centre = 256; ctx->pixel_size = 5e7; + ctx->reflectionlist = reflectionlist_new(); image_add(ctx->images, NULL, 512, 512, deg2rad(0.0), ctx); ctx->images->images[0].features = image_feature_list_new(); @@ -121,7 +122,7 @@ int main(int argc, char *argv[]) { 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); + ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real, ctx); // printf("RT: %i test features generated.\n", ctx->images->images[0].features->n_features); /* The "model" cell to be refined */ @@ -131,13 +132,13 @@ int main(int argc, char *argv[]) { ctx->cell->c.x = 0.1e9; ctx->cell->c.y = 0.1e9; ctx->cell->c.z = 5.3e9; ctx->cell_lattice = reflection_list_from_cell(ctx->cell); ctx->cell_lattice = reflection_list_from_cell(ctx->cell); - ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice); + ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice, ctx); 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); + ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice, ctx); fail = check_cell(ctx->cell, cell_real); @@ -157,6 +158,7 @@ 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) { }; +void displaywindow_message(DisplayWindow *dw, const char *text) { }; 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/refinetest3d1.c b/src/refinetest3d1.c index cd8cc84..4f7d706 100644 --- a/src/refinetest3d1.c +++ b/src/refinetest3d1.c @@ -3,7 +3,7 @@ * * Unit test for refinement procedure in 3D * - * (c) 2007-2008 Thomas White <taw27@cam.ac.uk> + * (c) 2007-2009 Thomas White <taw27@cam.ac.uk> * * dtr - Diffraction Tomography Reconstruction * @@ -27,10 +27,10 @@ #include "displaywindow.h" #include "mrc.h" -static int check_cell(Basis *cell, Basis *cell_real) { - +static int check_cell(Basis *cell, Basis *cell_real) +{ int fail; - + printf(" Refinement Actual\n"); printf("---------------------------\n"); printf("ax %+8f %+8f\n", cell->a.x/1e9, cell_real->a.x/1e9); @@ -42,129 +42,138 @@ static int check_cell(Basis *cell, Basis *cell_real) { printf("cx %+8f %+8f\n", cell->c.x/1e9, cell_real->c.x/1e9); printf("cy %+8f %+8f\n", cell->c.y/1e9, cell_real->c.y/1e9); printf("cz %+8f %+8f\n", cell->c.z/1e9, cell_real->c.z/1e9); - + fail = 0; if ( fabs(cell->a.x - cell_real->a.x) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: ax not determined correctly (got %8f, should be %8f)\n", - cell->a.x/1e9, cell_real->a.x/1e9); + fprintf(stderr, "refinetest3d: ax not determined correctly " + "(got %8f, should be %8f)\n", + cell->a.x/1e9, cell_real->a.x/1e9); fail = 1; } else { printf("ax is OK.\n"); } if ( fabs(cell->a.y - cell_real->a.y) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: ay not determined correctly (got %8f, should be %8f)\n", - cell->a.y/1e9, cell_real->a.y/1e9); + fprintf(stderr, "refinetest3d: ay not determined correctly " + "(got %8f, should be %8f)\n", + cell->a.y/1e9, cell_real->a.y/1e9); fail = 1; } else { printf("ay is OK.\n"); } if ( fabs(cell->a.z - cell_real->a.z) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: az not determined correctly (got %8f, should be %8f)\n", - cell->a.z/1e9, cell_real->a.z/1e9); + fprintf(stderr, "refinetest3d: az not determined correctly " + "(got %8f, should be %8f)\n", + cell->a.z/1e9, cell_real->a.z/1e9); fail = 1; } else { printf("az is OK.\n"); } if ( fabs(cell->b.x - cell_real->b.x) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: bx not determined correctly (got %8f, should be %8f)\n", - cell->b.x/1e9, cell_real->b.x/1e9); + fprintf(stderr, "refinetest3d: bx not determined correctly " + "(got %8f, should be %8f)\n", + cell->b.x/1e9, cell_real->b.x/1e9); fail = 1; } else { printf("bx is OK.\n"); } if ( fabs(cell->b.y - cell_real->b.y) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: by not determined correctly (got %8f, should be %8f)\n", - cell->b.y/1e9, cell_real->b.y/1e9); + fprintf(stderr, "refinetest3d: by not determined correctly " + "(got %8f, should be %8f)\n", + cell->b.y/1e9, cell_real->b.y/1e9); fail = 1; } else { printf("by is OK.\n"); } if ( fabs(cell->c.x - cell_real->c.x) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: cx not determined correctly (got %8f, should be %8f)\n", - cell->c.x/1e9, cell_real->c.x/1e9); + fprintf(stderr, "refinetest3d: cx not determined correctly " + "(got %8f, should be %8f)\n", + cell->c.x/1e9, cell_real->c.x/1e9); fail = 1; } else { printf("cx is OK.\n"); } if ( fabs(cell->c.z - cell_real->c.z) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: cz not determined correctly (got %8f, should be %8f)\n", - cell->c.z/1e9, cell_real->c.z/1e9); + fprintf(stderr, "refinetest3d: cz not determined correctly " + "(got %8f, should be %8f)\n", + cell->c.z/1e9, cell_real->c.z/1e9); fail = 1; } else { printf("cz is OK.\n"); } - - printf("bz and cy not checked in this test.\n"); - - return fail; -} + printf("bz and cy not checked in this test.\n"); -int main(int argc, char *argv[]) { + return fail; +} +int main(int argc, char *argv[]) +{ ControlContext *ctx; ReflectionList *reflections_real; Basis *cell_real; int fail; - + ctx = control_ctx_new(); - ctx->omega = deg2rad(0.0); + ctx->omega = deg2rad(45.0); ctx->lambda = lambda(300.0e3); /* 300 keV */ ctx->fmode = FORMULATION_PIXELSIZE; ctx->x_centre = 256; ctx->y_centre = 256; ctx->pixel_size = 5e7; + ctx->reflectionlist = reflectionlist_new(); image_add(ctx->images, NULL, 512, 512, deg2rad(00.0), ctx); ctx->images->images[0].features = image_feature_list_new(); image_add(ctx->images, NULL, 512, 512, deg2rad(90.0), ctx); ctx->images->images[1].features = image_feature_list_new(); - //ctx->omega = deg2rad(90.0); - //image_add(ctx->images, NULL, 512, 512, deg2rad(90.0), ctx); - //ctx->images->images[2].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; + 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); - ctx->images->images[1].features = reproject_get_reflections(&ctx->images->images[1], reflections_real); - //ctx->images->images[2].features = reproject_get_reflections(&ctx->images->images[2], reflections_real); - + ctx->images->images[0].features = + reproject_get_reflections(&ctx->images->images[0], + reflections_real, ctx); + ctx->images->images[1].features = + reproject_get_reflections(&ctx->images->images[1], + reflections_real, ctx); + /* The "model" cell to be refined */ ctx->cell = malloc(sizeof(Basis)); - ctx->cell->a.x = 5.0e9; ctx->cell->a.y = 0.1e9; ctx->cell->a.z = 0.1e9; - ctx->cell->b.x = 0.1e9; ctx->cell->b.y = 5.0e9; ctx->cell->b.z = 0.1e9; - ctx->cell->c.x = 0.1e9; ctx->cell->c.y = 0.1e9; ctx->cell->c.z = 5.0e9; + ctx->cell->a.x = 5.2e9; ctx->cell->a.y = 0.1e9; ctx->cell->a.z = 0.1e9; + ctx->cell->b.x = 0.2e9; ctx->cell->b.y = 4.8e9; ctx->cell->b.z = 0.1e9; + ctx->cell->c.x = 0.1e9; ctx->cell->c.y = 0.1e9; ctx->cell->c.z = 5.3e9; 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]); - ctx->images->images[1].rflist = reproject_get_reflections(&ctx->images->images[1], ctx->cell_lattice); - reproject_partner_features(ctx->images->images[1].rflist, &ctx->images->images[1]); - //ctx->images->images[2].rflist = reproject_get_reflections(&ctx->images->images[2], ctx->cell_lattice); - //reproject_partner_features(ctx->images->images[2].rflist, &ctx->images->images[2]); - + ctx->images->images[0].rflist = + reproject_get_reflections(&ctx->images->images[0], + ctx->cell_lattice, ctx); + reproject_partner_features(ctx->images->images[0].rflist, + &ctx->images->images[0]); + ctx->images->images[1].rflist = + reproject_get_reflections(&ctx->images->images[1], + ctx->cell_lattice, ctx); + reproject_partner_features(ctx->images->images[1].rflist, + &ctx->images->images[1]); + refine_do_cell(ctx); image_feature_list_free(ctx->images->images[0].rflist); image_feature_list_free(ctx->images->images[1].rflist); - //image_feature_list_free(ctx->images->images[2].rflist); - + fail = check_cell(ctx->cell, cell_real); - + free(ctx); - + if ( fail ) return 1; - + printf("\n3D refinement test OK.\n"); - + return 0; - } /* Dummy function stubs */ @@ -173,7 +182,8 @@ 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) { }; +void displaywindow_message(DisplayWindow *dw, const char *text) { }; 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) { }; - +void gtk_value_graph_set_data(GtkValueGraph *vg, double *data, unsigned int n) + { }; diff --git a/src/refinetest3d2.c b/src/refinetest3d2.c index 38ea72e..25a3f1f 100644 --- a/src/refinetest3d2.c +++ b/src/refinetest3d2.c @@ -3,7 +3,7 @@ * * Unit test for refinement procedure in 3D * - * (c) 2007-2008 Thomas White <taw27@cam.ac.uk> + * (c) 2007-2009 Thomas White <taw27@cam.ac.uk> * * dtr - Diffraction Tomography Reconstruction * @@ -28,9 +28,9 @@ #include "mrc.h" static int check_cell(Basis *cell, Basis *cell_real) { - + int fail; - + printf(" Refinement Actual\n"); printf("---------------------------\n"); printf("ax %+8f %+8f\n", cell->a.x/1e9, cell_real->a.x/1e9); @@ -42,75 +42,84 @@ static int check_cell(Basis *cell, Basis *cell_real) { printf("cx %+8f %+8f\n", cell->c.x/1e9, cell_real->c.x/1e9); printf("cy %+8f %+8f\n", cell->c.y/1e9, cell_real->c.y/1e9); printf("cz %+8f %+8f\n", cell->c.z/1e9, cell_real->c.z/1e9); - + fail = 0; if ( fabs(cell->a.x - cell_real->a.x) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: ax not determined correctly (got %8f, should be %8f)\n", - cell->a.x/1e9, cell_real->a.x/1e9); + fprintf(stderr, "refinetest3d: ax not determined correctly " + "(got %8f, should be %8f)\n", + cell->a.x/1e9, cell_real->a.x/1e9); fail = 1; } else { printf("ax is OK.\n"); } if ( fabs(cell->a.y - cell_real->a.y) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: ay not determined correctly (got %8f, should be %8f)\n", - cell->a.y/1e9, cell_real->a.y/1e9); + fprintf(stderr, "refinetest3d: ay not determined correctly " + "(got %8f, should be %8f)\n", + cell->a.y/1e9, cell_real->a.y/1e9); fail = 1; } else { printf("ay is OK.\n"); } if ( fabs(cell->a.z - cell_real->a.z) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: az not determined correctly (got %8f, should be %8f)\n", - cell->a.z/1e9, cell_real->a.z/1e9); + fprintf(stderr, "refinetest3d: az not determined correctly " + "(got %8f, should be %8f)\n", + cell->a.z/1e9, cell_real->a.z/1e9); fail = 1; } else { printf("az is OK.\n"); } if ( fabs(cell->b.x - cell_real->b.x) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: bx not determined correctly (got %8f, should be %8f)\n", - cell->b.x/1e9, cell_real->b.x/1e9); + fprintf(stderr, "refinetest3d: bx not determined correctly " + "(got %8f, should be %8f)\n", + cell->b.x/1e9, cell_real->b.x/1e9); fail = 1; } else { printf("bx is OK.\n"); } if ( fabs(cell->b.y - cell_real->b.y) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: by not determined correctly (got %8f, should be %8f)\n", - cell->b.y/1e9, cell_real->b.y/1e9); + fprintf(stderr, "refinetest3d: by not determined correctly " + "(got %8f, should be %8f)\n", + cell->b.y/1e9, cell_real->b.y/1e9); fail = 1; } else { printf("by is OK.\n"); } if ( fabs(cell->b.z - cell_real->b.z) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: bz not determined correctly (got %8f, should be %8f)\n", - cell->b.z/1e9, cell_real->b.z/1e9); + fprintf(stderr, "refinetest3d: bz not determined correctly " + "(got %8f, should be %8f)\n", + cell->b.z/1e9, cell_real->b.z/1e9); fail = 1; } else { printf("bz is OK.\n"); } if ( fabs(cell->c.x - cell_real->c.x) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: cx not determined correctly (got %8f, should be %8f)\n", - cell->c.x/1e9, cell_real->c.x/1e9); + fprintf(stderr, "refinetest3d: cx not determined correctly " + "(got %8f, should be %8f)\n", + cell->c.x/1e9, cell_real->c.x/1e9); fail = 1; } else { printf("cx is OK.\n"); } if ( fabs(cell->c.y - cell_real->c.y) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: cy not determined correctly (got %8f, should be %8f)\n", - cell->c.y/1e9, cell_real->c.y/1e9); + fprintf(stderr, "refinetest3d: cy not determined correctly " + "(got %8f, should be %8f)\n", + cell->c.y/1e9, cell_real->c.y/1e9); fail = 1; } else { printf("cx is OK.\n"); } if ( fabs(cell->c.z - cell_real->c.z) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: cz not determined correctly (got %8f, should be %8f)\n", - cell->c.z/1e9, cell_real->c.z/1e9); + fprintf(stderr, "refinetest3d: cz not determined correctly " + "(got %8f, should be %8f)\n", + cell->c.z/1e9, cell_real->c.z/1e9); fail = 1; } else { printf("cz is OK.\n"); } - + return fail; -} +} int main(int argc, char *argv[]) { @@ -118,7 +127,7 @@ int main(int argc, char *argv[]) { ReflectionList *reflections_real; Basis *cell_real; int fail; - + ctx = control_ctx_new(); ctx->omega = deg2rad(0.0); ctx->lambda = lambda(300.0e3); /* 300 keV */ @@ -126,6 +135,7 @@ int main(int argc, char *argv[]) { ctx->x_centre = 256; ctx->y_centre = 256; ctx->pixel_size = 5e7; + ctx->reflectionlist = reflectionlist_new(); image_add(ctx->images, NULL, 512, 512, deg2rad(00.0), ctx); ctx->images->images[0].features = image_feature_list_new(); image_add(ctx->images, NULL, 512, 512, deg2rad(90.0), ctx); @@ -133,50 +143,65 @@ int main(int argc, char *argv[]) { ctx->omega = deg2rad(90.0); image_add(ctx->images, NULL, 512, 512, deg2rad(90.0), ctx); ctx->images->images[2].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; + 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); - ctx->images->images[1].features = reproject_get_reflections(&ctx->images->images[1], reflections_real); - ctx->images->images[2].features = reproject_get_reflections(&ctx->images->images[2], reflections_real); - + ctx->images->images[0].features = + reproject_get_reflections(&ctx->images->images[0], + reflections_real, ctx); + ctx->images->images[1].features = + reproject_get_reflections(&ctx->images->images[1], + reflections_real, ctx); + ctx->images->images[2].features = + reproject_get_reflections(&ctx->images->images[2], + reflections_real, ctx); + /* The "model" cell to be refined */ ctx->cell = malloc(sizeof(Basis)); - ctx->cell->a.x = 5.0e9; ctx->cell->a.y = 0.1e9; ctx->cell->a.z = 0.1e9; - ctx->cell->b.x = 0.1e9; ctx->cell->b.y = 5.0e9; ctx->cell->b.z = 0.1e9; - ctx->cell->c.x = 0.1e9; ctx->cell->c.y = 0.1e9; ctx->cell->c.z = 5.0e9; + ctx->cell->a.x = 5.2e9; ctx->cell->a.y = 0.1e9; ctx->cell->a.z = 0.1e9; + ctx->cell->b.x = 0.2e9; ctx->cell->b.y = 4.8e9; ctx->cell->b.z = 0.1e9; + ctx->cell->c.x = 0.1e9; ctx->cell->c.y = 0.1e9; ctx->cell->c.z = 5.3e9; 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]); - ctx->images->images[1].rflist = reproject_get_reflections(&ctx->images->images[1], ctx->cell_lattice); - reproject_partner_features(ctx->images->images[1].rflist, &ctx->images->images[1]); - ctx->images->images[2].rflist = reproject_get_reflections(&ctx->images->images[2], ctx->cell_lattice); - reproject_partner_features(ctx->images->images[2].rflist, &ctx->images->images[2]); - + ctx->images->images[0].rflist = + reproject_get_reflections(&ctx->images->images[0], + ctx->cell_lattice, ctx); + reproject_partner_features(ctx->images->images[0].rflist, + &ctx->images->images[0]); + ctx->images->images[1].rflist = + reproject_get_reflections(&ctx->images->images[1], + ctx->cell_lattice, ctx); + reproject_partner_features(ctx->images->images[1].rflist, + &ctx->images->images[1]); + ctx->images->images[2].rflist = + reproject_get_reflections(&ctx->images->images[2], + ctx->cell_lattice, ctx); + reproject_partner_features(ctx->images->images[2].rflist, + &ctx->images->images[2]); + refine_do_cell(ctx); image_feature_list_free(ctx->images->images[0].rflist); image_feature_list_free(ctx->images->images[1].rflist); image_feature_list_free(ctx->images->images[2].rflist); - + fail = check_cell(ctx->cell, cell_real); - + free(ctx); - + if ( fail ) return 1; - + printf("\n3D refinement test OK.\n"); - + return 0; - + } /* Dummy function stubs */ @@ -185,7 +210,8 @@ 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) { }; +void displaywindow_message(DisplayWindow *dw, const char *text) { }; 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) { }; - +void gtk_value_graph_set_data(GtkValueGraph *vg, double *data, unsigned int n) + { }; diff --git a/src/refinetest3d3.c b/src/refinetest3d3.c index de8a792..1a1623d 100644 --- a/src/refinetest3d3.c +++ b/src/refinetest3d3.c @@ -3,7 +3,7 @@ * * Unit test for refinement procedure in 3D * - * (c) 2007-2008 Thomas White <taw27@cam.ac.uk> + * (c) 2007-2009 Thomas White <taw27@cam.ac.uk> * * dtr - Diffraction Tomography Reconstruction * @@ -28,9 +28,9 @@ #include "mrc.h" static int check_cell(Basis *cell, Basis *cell_real) { - + int fail; - + printf(" Refinement Actual\n"); printf("---------------------------\n"); printf("ax %+8f %+8f\n", cell->a.x/1e9, cell_real->a.x/1e9); @@ -42,75 +42,84 @@ static int check_cell(Basis *cell, Basis *cell_real) { printf("cx %+8f %+8f\n", cell->c.x/1e9, cell_real->c.x/1e9); printf("cy %+8f %+8f\n", cell->c.y/1e9, cell_real->c.y/1e9); printf("cz %+8f %+8f\n", cell->c.z/1e9, cell_real->c.z/1e9); - + fail = 0; if ( fabs(cell->a.x - cell_real->a.x) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: ax not determined correctly (got %8f, should be %8f)\n", - cell->a.x/1e9, cell_real->a.x/1e9); + fprintf(stderr, "refinetest3d: ax not determined correctly " + "(got %8f, should be %8f)\n", + cell->a.x/1e9, cell_real->a.x/1e9); fail = 1; } else { printf("ax is OK.\n"); } if ( fabs(cell->a.y - cell_real->a.y) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: ay not determined correctly (got %8f, should be %8f)\n", - cell->a.y/1e9, cell_real->a.y/1e9); + fprintf(stderr, "refinetest3d: ay not determined correctly " + "(got %8f, should be %8f)\n", + cell->a.y/1e9, cell_real->a.y/1e9); fail = 1; } else { printf("ay is OK.\n"); } if ( fabs(cell->a.z - cell_real->a.z) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: az not determined correctly (got %8f, should be %8f)\n", - cell->a.z/1e9, cell_real->a.z/1e9); + fprintf(stderr, "refinetest3d: az not determined correctly " + "(got %8f, should be %8f)\n", + cell->a.z/1e9, cell_real->a.z/1e9); fail = 1; } else { printf("az is OK.\n"); } if ( fabs(cell->b.x - cell_real->b.x) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: bx not determined correctly (got %8f, should be %8f)\n", - cell->b.x/1e9, cell_real->b.x/1e9); + fprintf(stderr, "refinetest3d: bx not determined correctly " + "(got %8f, should be %8f)\n", + cell->b.x/1e9, cell_real->b.x/1e9); fail = 1; } else { printf("bx is OK.\n"); } if ( fabs(cell->b.y - cell_real->b.y) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: by not determined correctly (got %8f, should be %8f)\n", - cell->b.y/1e9, cell_real->b.y/1e9); + fprintf(stderr, "refinetest3d: by not determined correctly " + "(got %8f, should be %8f)\n", + cell->b.y/1e9, cell_real->b.y/1e9); fail = 1; } else { printf("by is OK.\n"); } if ( fabs(cell->b.z - cell_real->b.z) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: bz not determined correctly (got %8f, should be %8f)\n", - cell->b.z/1e9, cell_real->b.z/1e9); + fprintf(stderr, "refinetest3d: bz not determined correctly " + "(got %8f, should be %8f)\n", + cell->b.z/1e9, cell_real->b.z/1e9); fail = 1; } else { printf("bz is OK.\n"); } if ( fabs(cell->c.x - cell_real->c.x) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: cx not determined correctly (got %8f, should be %8f)\n", - cell->c.x/1e9, cell_real->c.x/1e9); + fprintf(stderr, "refinetest3d: cx not determined correctly " + "(got %8f, should be %8f)\n", + cell->c.x/1e9, cell_real->c.x/1e9); fail = 1; } else { printf("cx is OK.\n"); } if ( fabs(cell->c.y - cell_real->c.y) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: cy not determined correctly (got %8f, should be %8f)\n", - cell->c.y/1e9, cell_real->c.y/1e9); + fprintf(stderr, "refinetest3d: cy not determined correctly " + "(got %8f, should be %8f)\n", + cell->c.y/1e9, cell_real->c.y/1e9); fail = 1; } else { printf("cx is OK.\n"); } if ( fabs(cell->c.z - cell_real->c.z) > 0.01e9 ) { - fprintf(stderr, "refinetest3d: cz not determined correctly (got %8f, should be %8f)\n", - cell->c.z/1e9, cell_real->c.z/1e9); + fprintf(stderr, "refinetest3d: cz not determined correctly " + "(got %8f, should be %8f)\n", + cell->c.z/1e9, cell_real->c.z/1e9); fail = 1; } else { printf("cz is OK.\n"); } - + return fail; -} +} int main(int argc, char *argv[]) { @@ -119,7 +128,7 @@ int main(int argc, char *argv[]) { Basis *cell_real; int fail; int i; - + ctx = control_ctx_new(); ctx->omega = deg2rad(45.0); ctx->lambda = lambda(300.0e3); /* 300 keV */ @@ -127,52 +136,58 @@ int main(int argc, char *argv[]) { ctx->x_centre = 256; ctx->y_centre = 256; ctx->pixel_size = 5e7; + ctx->reflectionlist = reflectionlist_new(); for ( i=0; i<=90; i++ ) { image_add(ctx->images, NULL, 512, 512, deg2rad(i), ctx); ctx->images->images[i].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; + 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); for ( i=0; i<=90; i++ ) { - ctx->images->images[i].features = reproject_get_reflections(&ctx->images->images[i], reflections_real); + ctx->images->images[i].features = + reproject_get_reflections(&ctx->images->images[i], + reflections_real, ctx); } - + /* The "model" cell to be refined */ ctx->cell = malloc(sizeof(Basis)); - ctx->cell->a.x = 5.0e9; ctx->cell->a.y = 0.1e9; ctx->cell->a.z = 0.1e9; - ctx->cell->b.x = 0.1e9; ctx->cell->b.y = 5.0e9; ctx->cell->b.z = 0.1e9; - ctx->cell->c.x = 0.1e9; ctx->cell->c.y = 0.1e9; ctx->cell->c.z = 5.0e9; + ctx->cell->a.x = 5.2e9; ctx->cell->a.y = 0.1e9; ctx->cell->a.z = 0.1e9; + ctx->cell->b.x = 0.2e9; ctx->cell->b.y = 4.8e9; ctx->cell->b.z = 0.1e9; + ctx->cell->c.x = 0.1e9; ctx->cell->c.y = 0.1e9; ctx->cell->c.z = 5.3e9; ctx->cell_lattice = reflection_list_from_cell(ctx->cell); for ( i=0; i<ctx->images->n_images; i++ ) { - ctx->images->images[i].rflist = reproject_get_reflections(&ctx->images->images[i], ctx->cell_lattice); - reproject_partner_features(ctx->images->images[i].rflist, &ctx->images->images[i]); + ctx->images->images[i].rflist = + reproject_get_reflections(&ctx->images->images[i], + ctx->cell_lattice, ctx); + reproject_partner_features(ctx->images->images[i].rflist, + &ctx->images->images[i]); } - + refine_do_cell(ctx); for ( i=0; i<ctx->images->n_images; i++ ) { image_feature_list_free(ctx->images->images[i].rflist); } - + fail = check_cell(ctx->cell, cell_real); - + free(ctx); - + if ( fail ) return 1; - + printf("\n3D refinement test OK.\n"); - + return 0; - + } /* Dummy function stubs */ @@ -181,7 +196,8 @@ 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) { }; +void displaywindow_message(DisplayWindow *dw, const char *text) { }; 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) { }; - +void gtk_value_graph_set_data(GtkValueGraph *vg, double *data, unsigned int n) + { }; diff --git a/src/reflections.h b/src/reflections.h index 8fafe4b..e51c530 100644 --- a/src/reflections.h +++ b/src/reflections.h @@ -22,7 +22,8 @@ typedef enum { REFLECTION_MARKER, REFLECTION_VECTOR_MARKER_1, REFLECTION_VECTOR_MARKER_2, - REFLECTION_VECTOR_MARKER_3 + REFLECTION_VECTOR_MARKER_3, + REFLECTION_VECTOR_MARKER_4 } ReflectionType; typedef struct reflection_struct { diff --git a/src/reproject.h b/src/reproject.h index 93f010a..37eb84d 100644 --- a/src/reproject.h +++ b/src/reproject.h @@ -19,7 +19,8 @@ #include "control.h" #include "image.h" -extern ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist); +extern ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist, + ControlContext *ctx); extern void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image); extern void reproject_cell_to_lattice(ControlContext *ctx); extern void reproject_lattice_changed(ControlContext *ctx); |