From d2a3b45e9624801b55f5d4158ecf885c2442b8c9 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 1 Apr 2009 17:50:03 +0100 Subject: Refinement tweaks and tidy-up --- src/refine.c | 278 +++++++++++++++++++++++++++++----------------------- src/refinetest3d1.c | 124 ++++++++++++----------- src/refinetest3d2.c | 128 ++++++++++++++---------- src/refinetest3d3.c | 108 +++++++++++--------- 4 files changed, 360 insertions(+), 278 deletions(-) diff --git a/src/refine.c b/src/refine.c index 19f9ea3..a642096 100644 --- a/src/refine.c +++ b/src/refine.c @@ -3,7 +3,7 @@ * * Refine the reconstruction * - * (c) 2007-2008 Thomas White + * (c) 2007-2009 Thomas White * * dtr - Diffraction Tomography Reconstruction * @@ -41,34 +41,38 @@ #define DELTA (0.0001e9) typedef struct { - signed int h; signed int k; signed int l; /* Indicies of reflection */ - double dx; double dy; /* Deviation vector for this reflection in nm^-1 */ - double xx; double xy; double xz; /* x,y,z components of "image x" */ - double yx; double yy; double yz; /* x,y,z components of "image y" */ + /* Indicies of reflection */ + signed int h; signed int k; signed int l; + /* 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 - int image; /* Image number this dev came from, for debugging purposes */ + /* Image number this dev came from, for debugging purposes */ + int image; double rdx; double rdy; #endif /* REFINE_DEBUG_MORE */ - -} Deviation; -static void refine_printf(DisplayWindow *dw, char *fmt, ...) { +} Deviation; +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) { - +void refine_do_sequence(ControlContext *ctx) +{ double omega_offs; int idx; double *fit_vals; @@ -76,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; iimages->n_images; i++ ) { ctx->images->images[i].omega += omega_offs; } reproject_lattice_changed(ctx); - + fit = refine_do_cell(ctx); - refine_printf(ctx->dw, " omega_offs=%f deg, fit=%f nm^-1\n", + 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; iimages->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"); @@ -118,7 +124,7 @@ 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 */ refine_printf(ctx->dw, "Best omega offset = %f deg (%f nm^-1)\n", rad2deg(omega_offs_best), fit_best/DISPFACTOR); @@ -128,50 +134,53 @@ void refine_do_sequence(ControlContext *ctx) { refine_do_cell(ctx); reproject_lattice_changed(ctx); mapping_adjust_axis(ctx, omega_offs_best); - } -static double refine_mean_dev(Deviation *d, int nf, Basis *devcell, int disp) { - +static double refine_mean_dev(Deviation *d, int nf, Basis *devcell, int disp) +{ double fom = 0.0; int f; - + for ( f=0; fa.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; - + 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; + /* 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; - + dx -= d[f].dx; dy -= d[f].dy; - + #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); + 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); - + } - + 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; @@ -192,20 +201,24 @@ static void refine_cell_delta(Basis *devcell_try, int comp) { 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"); + default : fprintf(stderr, + "refine_cell_delta: argument out of range\n"); } - + } 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); + 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); + 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); - + cell.c.x/DISPFACTOR, cell.c.y/DISPFACTOR, + cell.c.z/DISPFACTOR); + } static Basis refine_test_confidence(Deviation *d, int nf) { @@ -213,58 +226,58 @@ 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; } @@ -281,66 +294,72 @@ double refine_do_cell(ControlContext *ctx) { 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; iimages->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); + ctx->images->images[i].rflist = + reproject_get_reflections(&ctx->images->images[i], + ctx->cell_lattice, ctx); } - + for ( j=0; jimages->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++; } - + refine_printf(ctx->dw, "%i features from image %i\n", nf, i); - + } - refine_printf(ctx->dw, "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; iimages->n_images; i++ ) { - + ImageRecord *image; int j; - + image = &ctx->images->images[i]; - + for ( j=0; jimages->images[i].rflist->n_features; j++ ) { - + ImageFeature *rf; double dix, diy, dx, dy; 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; - 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); - + 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; @@ -349,94 +368,114 @@ double refine_do_cell(ControlContext *ctx) { rf->partner->x = old_x; rf->partner->y = old_y; double mod = sqrt(dx*dx + dy*dy)/DISPFACTOR; - refine_printf(ctx->dw, "= %+10.5f %+10.5f nm^-1 (length %7.5f nm^1)\n", - dx/DISPFACTOR, dy/DISPFACTOR, mod); - + 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); - + 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 ); - + /* 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); - + 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; - + #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, "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); + 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); + 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); + 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++ ) { - + memcpy(devcell_try, devcell, sizeof(Basis)); refine_cell_delta(devcell_try, comp1); refine_cell_delta(devcell_try, comp2); - mean_dev_try = refine_mean_dev(d, nf, devcell_try, 0); - + mean_dev_try = refine_mean_dev(d, nf, + devcell_try, 0); + /* Improvement greater than the tolerance? */ - if ( mean_dev_try < mean_dev ) { + if ( mean_dev_try < mean_dev-0.001 ) { mean_dev = mean_dev_try; - memcpy(devcell_try_best, devcell_try, sizeof(Basis)); + 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); + refine_printf(ctx->dw, + "After %5i iterations: mean dev = %13.8f nm^1\n", + it, mean_dev/DISPFACTOR); } - + } 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 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 */ conf_threshold = 2.5; if ( conf.a.x > conf_threshold ) ctx->cell->a.x += devcell->a.x; @@ -448,12 +487,11 @@ double refine_do_cell(ControlContext *ctx) { 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 mean_dev; - -} +} diff --git a/src/refinetest3d1.c b/src/refinetest3d1.c index 6dfc74a..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 + * (c) 2007-2009 Thomas White * * 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,73 +42,79 @@ 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; @@ -119,53 +125,55 @@ int main(int argc, char *argv[]) { 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); - 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->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, 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); - //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 */ @@ -177,5 +185,5 @@ 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 cffe7f2..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 + * (c) 2007-2009 Thomas White * * 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 */ @@ -134,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); - 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); - + 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, 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]); - + 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 */ @@ -189,5 +213,5 @@ 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 65705a8..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 + * (c) 2007-2009 Thomas White * * 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 */ @@ -132,50 +141,53 @@ int main(int argc, char *argv[]) { 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); + 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; iimages->n_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]); + 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; iimages->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 */ @@ -187,5 +199,5 @@ 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) + { }; -- cgit v1.2.3