diff options
Diffstat (limited to 'src/refine.c')
-rw-r--r-- | src/refine.c | 100 |
1 files changed, 75 insertions, 25 deletions
diff --git a/src/refine.c b/src/refine.c index a5ae3d5..7838e56 100644 --- a/src/refine.c +++ b/src/refine.c @@ -27,6 +27,7 @@ #include "reproject.h" #include "control.h" #include "mapping.h" +#include "imagedisplay.h" /* Return the root sum squared deviation distance for all the "reprojectable" features in an image */ static double refine_image_deviation(ImageRecord *image, ReflectionList *cell_lattice) { @@ -63,16 +64,36 @@ typedef enum { INDEX_C = 1<<2 } RefinementIndex; +static const char *refine_decode(RefinementIndex i) { + + switch ( i ) { + + case INDEX_A : return "a--"; + case INDEX_B : return "-b-"; + case INDEX_C : return "--c"; + case INDEX_A | INDEX_B : return "ab-"; + case INDEX_A | INDEX_C : return "a-c"; + case INDEX_B | INDEX_C : return "-bc"; + case INDEX_A | INDEX_B | INDEX_C : return "abc"; + + default : return "---"; + + } + +} + /* Use the IPR algorithm to make "cell" fit the given image */ -static void refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *cell_lattice) { +static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *cell_lattice) { ImageFeatureList *rflist; ImageFeatureList *flist; + ImageFeature *ret; int i; Basis cd; /* Cell delta */ int n_a = 0; int n_b = 0; int n_c = 0; + int done = FALSE; cd.a.x = 0.0; cd.a.y = 0.0; cd.a.z = 0.0; cd.b.x = 0.0; cd.b.y = 0.0; cd.b.z = 0.0; @@ -131,54 +152,71 @@ static void refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *ce dk = ((a13*a32-a12*a33)*dx + (a11*a33-a13*a31)*dy + (a12*a31-a11*a32)*dz) / det; dl = ((a12*a23-a13*a22)*dx + (a13*a21-a11*a23)*dy + (a11*a22-a12*a21)*dz) / det; printf("dev(hkl) = %f %f %f\n", dh, dk, dl); - + delta = 0; - if ( fabs(dh) < 0.001 ) delta |= INDEX_A; - if ( fabs(dk) < 0.001 ) delta |= INDEX_B; - if ( fabs(dl) < 0.001 ) delta |= INDEX_C; - /* Typically, 'delta' will have all three bits set */ + if ( fabs(dh)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_A; + if ( fabs(dk)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_B; + if ( fabs(dl)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_C; shared = index & delta; + printf(" index: %s\n delta: %s\nshared: %s\n", refine_decode(index), refine_decode(delta), refine_decode(shared)); + if ( shared == 0 ) { /* No indices left - 'pure shear' (delta is perpendicular (in the abc basis) to index) */ shared = index; + printf("Pure shear.\n"); } if ( shared & INDEX_A ) { - double w = abs(h) / (abs(h)+abs(k)+abs(l)); - cd.a.x += w*dx / h; - cd.a.y += w*dy / h; - cd.a.z += w*dz / h; + double w = (double)abs(h) / (abs(h)+abs(k)+abs(l)); + printf("w(a) = %f\n", w); + cd.a.x += w*dx / (double)h; + cd.a.y += w*dy / (double)h; + cd.a.z += w*dz / (double)h; n_a++; } if ( shared & INDEX_B ) { - double w = abs(k) / (abs(h)+abs(k)+abs(l)); - cd.b.x += w*dx / k; - cd.b.y += w*dy / k; - cd.b.z += w*dz / k; + double w = (double)abs(k) / (abs(h)+abs(k)+abs(l)); + printf("w(b) = %f\n", w); + cd.b.x += w*dx / (double)k; + cd.b.y += w*dy / (double)k; + cd.b.z += w*dz / (double)k; n_b++; } if ( shared & INDEX_C ) { - double w = abs(l) / (abs(h)+abs(k)+abs(l)); - cd.c.x += w*dx / l; - cd.c.y += w*dy / l; - cd.c.z += w*dz / l; + double w = (double)abs(l) / (abs(h)+abs(k)+abs(l)); + printf("w(c) = %f\n", w); + cd.c.x += w*dx / (double)l; + cd.c.y += w*dy / (double)l; + cd.c.z += w*dz / (double)l; n_c++; } - + + printf("Distortion along x: %+8e = %+8e\n", h*cd.a.x + k*cd.b.x + l*cd.c.x, dx); + printf("Distortion along y: %+8e = %+8e\n", h*cd.a.y + k*cd.b.y + l*cd.c.y, dy); + printf("Distortion along z: %+8e = %+8e\n", h*cd.a.z + k*cd.b.z + l*cd.c.z, dz); + + done = TRUE; + + break; + } - image_feature_list_free(rflist); + if ( !done ) { + printf("RF: No partners found.\n"); + return NULL; + } + printf("n_a=%i, n_b=%i, n_c=%i\n", n_a, n_b, n_c); if ( n_a ) { - cd.a.x /= n_a; cd.a.y /= n_a; cd.a.z /= n_a; + cd.a.x /= (double)n_a; cd.a.y /= (double)n_a; cd.a.z /= (double)n_a; } if ( n_b ) { - cd.b.x /= n_b; cd.b.y /= n_b; cd.b.z /= n_b; + cd.b.x /= (double)n_b; cd.b.y /= (double)n_b; cd.b.z /= (double)n_b; } if ( n_c ) { - cd.c.x /= n_c; cd.c.y /= n_c; cd.c.z /= n_c; + cd.c.x /= (double)n_c; cd.c.y /= (double)n_c; cd.c.z /= (double)n_c; } printf("Total distortion(a) = %+8e %+8e %+8e\n", cd.a.x, cd.a.y, cd.a.z); @@ -189,6 +227,11 @@ static void refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *ce cell->b.x += cd.b.x; cell->b.y += cd.b.y; cell->b.z += cd.b.z; cell->c.x += cd.c.x; cell->c.y += cd.c.y; cell->c.z += cd.c.z; + ret = rflist->features[i].partner; + image_feature_list_free(rflist); + + return ret; + } /* Display a graph of root sum squared deviation distance against some other parameter */ @@ -231,8 +274,15 @@ static gint refine_graph(GtkWidget *step_button, ControlContext *ctx) { static gint refine_step(GtkWidget *step_button, ControlContext *ctx) { if ( (ctx->reproject_id) && (ctx->cell_lattice) ) { - refine_fit_image(ctx->cell, &ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice); - reproject_lattice_changed(ctx); + + ImageFeature *fitted; + + fitted = refine_fit_image(ctx->cell, &ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice); + if ( fitted ) { + reproject_lattice_changed(ctx); + imagedisplay_add_mark(ctx->reproject_id, fitted->x,fitted->y, IMAGEDISPLAY_MARK_CIRCLE_3); + } + } else { displaywindow_error("Please first open the reprojection window and select the image to fit", ctx->dw); } |