From 4b6725f92bd05657134476056c88c3e224c10101 Mon Sep 17 00:00:00 2001 From: taw27 Date: Fri, 2 Nov 2007 19:06:04 +0000 Subject: More work on refinement Store absolute image coordinates in cache (not relative to centre) git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@190 bf6ca9ba-c028-0410-8290-897cf20841d1 --- src/basis.c | 2 +- src/basis.h | 2 + src/imagedisplay.c | 9 ++++- src/imagedisplay.h | 2 + src/itrans-stat.c | 2 +- src/itrans-threshold.c | 4 +- src/itrans-zaefferer.c | 3 +- src/mapping.c | 6 +-- src/refine.c | 100 ++++++++++++++++++++++++++++++++++++------------- src/reproject.c | 22 +++++------ 10 files changed, 104 insertions(+), 48 deletions(-) diff --git a/src/basis.c b/src/basis.c index 841ccc6..08e3edb 100644 --- a/src/basis.c +++ b/src/basis.c @@ -18,7 +18,7 @@ #include "reflections.h" #include "basis.h" -static double basis_efom(ReflectionList *reflectionlist, Basis *basis) { +double basis_efom(ReflectionList *reflectionlist, Basis *basis) { int n_indexed, n_counted; Reflection *cur; diff --git a/src/basis.h b/src/basis.h index 9fb1cda..5803cab 100644 --- a/src/basis.h +++ b/src/basis.h @@ -32,5 +32,7 @@ typedef struct basis_struct { Vector c; } Basis; +extern double basis_efom(struct reflectionlist_struct *reflectionlist, Basis *basis); + #endif /* BASIS_H */ diff --git a/src/imagedisplay.c b/src/imagedisplay.c index 3877da8..ceb36ff 100644 --- a/src/imagedisplay.c +++ b/src/imagedisplay.c @@ -210,11 +210,14 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even switch ( cur->type ) { case IMAGEDISPLAY_MARK_CIRCLE_1 : gc = imagedisplay->gc_marks_1; break; case IMAGEDISPLAY_MARK_CIRCLE_2 : gc = imagedisplay->gc_marks_2; break; + case IMAGEDISPLAY_MARK_CIRCLE_3 : gc = imagedisplay->gc_marks_3; break; case IMAGEDISPLAY_MARK_LINE_1 : gc = imagedisplay->gc_marks_1; break; + case IMAGEDISPLAY_MARK_LINE_2 : gc = imagedisplay->gc_marks_2; break; default : gc = imagedisplay->gc_marks_1; break; } - if ( (cur->type == IMAGEDISPLAY_MARK_CIRCLE_1) || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_2) ) { + if ( (cur->type == IMAGEDISPLAY_MARK_CIRCLE_1) || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_2) + || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_3) ) { gdk_draw_arc(drawingarea->window, gc, FALSE, xoffs + cur->x*scale - 5, @@ -258,6 +261,10 @@ static gint imagedisplay_realize(GtkWidget *widget, ImageDisplay *imagedisplay) gdk_color_parse("#00dd00", &colour); gdk_gc_set_rgb_fg_color(imagedisplay->gc_marks_2, &colour); + imagedisplay->gc_marks_3 = gdk_gc_new(imagedisplay->drawingarea->window); + gdk_color_parse("#00ddff", &colour); + gdk_gc_set_rgb_fg_color(imagedisplay->gc_marks_3, &colour); + imagedisplay->realised = TRUE; return 0; diff --git a/src/imagedisplay.h b/src/imagedisplay.h index 3b7192a..d404bcb 100644 --- a/src/imagedisplay.h +++ b/src/imagedisplay.h @@ -31,6 +31,7 @@ typedef enum { typedef enum { IMAGEDISPLAY_MARK_CIRCLE_1, IMAGEDISPLAY_MARK_CIRCLE_2, + IMAGEDISPLAY_MARK_CIRCLE_3, IMAGEDISPLAY_MARK_LINE_1, IMAGEDISPLAY_MARK_LINE_2 } ImageDisplayMarkType; @@ -66,6 +67,7 @@ typedef struct imagedisplay_struct { GdkGC *gc_tiltaxis; GdkGC *gc_marks_1; GdkGC *gc_marks_2; + GdkGC *gc_marks_3; gboolean realised; unsigned int drawingarea_width; diff --git a/src/itrans-stat.c b/src/itrans-stat.c index e08bd26..67ed96b 100644 --- a/src/itrans-stat.c +++ b/src/itrans-stat.c @@ -519,7 +519,7 @@ ImageFeatureList *itrans_peaksearch_stat(ImageRecord *imagerecord) { px = gsl_matrix_get(p,0,i); py = gsl_matrix_get(p,1,i); - image_add_feature(flist, (px-imagerecord->x_centre), (py-imagerecord->y_centre), imagerecord, 1.0); + image_add_feature(flist, px, py, imagerecord, 1.0); } diff --git a/src/itrans-threshold.c b/src/itrans-threshold.c index 17a2cbe..87a461f 100644 --- a/src/itrans-threshold.c +++ b/src/itrans-threshold.c @@ -45,7 +45,7 @@ ImageFeatureList *itrans_peaksearch_threshold(ImageRecord *imagerecord, ControlC assert(y=0); assert(y>=0); - image_add_feature(flist, (x-imagerecord->x_centre), (y-imagerecord->y_centre), imagerecord, image[x + width*y]); + image_add_feature(flist, x, y, imagerecord, image[x + width*y]); } } } @@ -100,7 +100,7 @@ ImageFeatureList *itrans_peaksearch_adaptive_threshold(ImageRecord *imagerecord, assert(max_y=0); assert(max_y>=0); - image_add_feature(flist, (x-imagerecord->x_centre), (y-imagerecord->y_centre), imagerecord, image[x + width*y]); + image_add_feature(flist, x, y, imagerecord, image[x + width*y]); /* Remove it and its surroundings */ for ( y=((max_y-10>0)?max_y-10:0); y<((max_y+10)=0); assert(mask_y>=0); - image_add_feature(flist, (mask_x-imagerecord->x_centre), (mask_y-imagerecord->y_centre), - imagerecord, image[mask_x + width*mask_y]); + image_add_feature(flist, mask_x, mask_y, imagerecord, image[mask_x + width*mask_y]); } } diff --git a/src/mapping.c b/src/mapping.c index 668bfc4..afe36c4 100644 --- a/src/mapping.c +++ b/src/mapping.c @@ -38,7 +38,7 @@ int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *d double nx, ny, nz; imagerecord = refl->parent; - x = refl->x; y = refl->y; + x = refl->x - imagerecord->x_centre; y = refl->y - imagerecord->y_centre; tilt = 2*M_PI*(imagerecord->tilt/360); /* Convert to Radians */ omega = 2*M_PI*(imagerecord->omega/360); /* Likewise */ k = 1/imagerecord->lambda; @@ -95,7 +95,7 @@ int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *d static void mapping_map_features(ControlContext *ctx) { int i; - + /* Create reflection list for measured reflections */ if ( ctx->reflectionlist ) reflectionlist_free(ctx->reflectionlist); ctx->reflectionlist = reflectionlist_new(); @@ -114,7 +114,7 @@ static void mapping_map_features(ControlContext *ctx) { reflection_add(ctx->reflectionlist, nx, ny, nz, ctx->images->images[i].features->features[j].intensity, REFLECTION_NORMAL); } - + } } 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); } diff --git a/src/reproject.c b/src/reproject.c index c0dfc73..2ff7e6c 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -85,7 +85,6 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * double theta; double x, y; double rx, ry, rz; - double xtest, ytest; /* Determine the intersection point */ xi = xl + s*nx; yi = yl + s*ny; zi = zl + s*nz; @@ -154,10 +153,11 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * return NULL; } + x += image->x_centre; + y += image->y_centre; + /* Sanity check */ - xtest = x + image->x_centre; - ytest = y + image->y_centre; - if ( (xtest>=0) && (xtestwidth) && (ytest>=0) && (ytestheight) ) { + if ( (x>=0) && (xwidth) && (y>=0) && (yheight) ) { /* Record the reflection */ image_add_feature_index(flist, x, y, image, reflection->intensity, @@ -171,7 +171,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * } reflection = reflection->next; - + } while ( reflection ); return flist; @@ -229,29 +229,25 @@ static void reproject_mark_peaks(ControlContext *ctx) { ImageFeatureList *rflist; ImageFeatureList *flist; size_t j; - double xc, yc; - - xc = ctx->images->images[ctx->reproject_cur_image].x_centre; - yc = ctx->images->images[ctx->reproject_cur_image].y_centre; /* Draw the reprojected peaks */ rflist = reproject_get_reflections(&ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice); for ( j=0; jn_features; j++ ) { - imagedisplay_add_mark(ctx->reproject_id, xc+rflist->features[j].x, yc+rflist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_1); + imagedisplay_add_mark(ctx->reproject_id, rflist->features[j].x, rflist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_1); } /* Now draw the original measured peaks */ flist = ctx->images->images[ctx->reproject_cur_image].features; for ( j=0; jn_features; j++ ) { - imagedisplay_add_mark(ctx->reproject_id, xc+flist->features[j].x, yc+flist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_2); + imagedisplay_add_mark(ctx->reproject_id, flist->features[j].x, flist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_2); } /* Now connect partners */ reproject_partner_features(rflist, &ctx->images->images[ctx->reproject_cur_image]); for ( j=0; jn_features; j++ ) { if ( rflist->features[j].partner ) { - imagedisplay_add_line(ctx->reproject_id, xc+rflist->features[j].x, yc+rflist->features[j].y, - xc+rflist->features[j].partner->x, yc+rflist->features[j].partner->y, + imagedisplay_add_line(ctx->reproject_id, rflist->features[j].x, rflist->features[j].y, + rflist->features[j].partner->x, rflist->features[j].partner->y, IMAGEDISPLAY_MARK_LINE_1); } } -- cgit v1.2.3