aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/displaywindow.c13
-rw-r--r--src/displaywindow.h1
-rw-r--r--src/intensities.c2
-rw-r--r--src/refine.c634
-rw-r--r--src/refinetest2d.c8
-rw-r--r--src/refinetest3d1.c126
-rw-r--r--src/refinetest3d2.c130
-rw-r--r--src/refinetest3d3.c108
-rw-r--r--src/reflections.h3
-rw-r--r--src/reproject.h3
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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-28 22:45:22 +0000