From 6a9b94d74b0d96f901d7a6baa885c306151ae0f4 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 8 Apr 2009 20:12:27 +0100 Subject: More de-clobbering of reprojection --- src/reproject.c | 206 ++++++++++++++++++++++++++++++++++++-------------------- 1 file changed, 132 insertions(+), 74 deletions(-) diff --git a/src/reproject.c b/src/reproject.c index df85f4f..e0199af 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -18,6 +18,7 @@ #include "imagedisplay.h" #include "displaywindow.h" #include "image.h" +#include "mapping.h" /* Attempt to find partners from the feature list of "image" for each feature in "flist". */ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) { @@ -37,6 +38,8 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) { if ( (d <= 20.0) && partner ) { rflist->features[i].partner = partner; rflist->features[i].partner_d = d; + image->features->features[idx].partner = + &rflist->features[i]; } else { rflist->features[i].partner = NULL; } @@ -45,16 +48,19 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) { } -ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist) { - +ImageFeatureList *reproject_get_reflections(ImageRecord *image, + ReflectionList *reflectionlist) +{ ImageFeatureList *flist; Reflection *reflection; - double smax = 0.4e9; + double smax = 0.2e9; double tilt, omega, k; - double nx, ny, nz; /* "normal" vector */ - double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda */ - double ux, uy, uz, uxt, uyt, uzt; /* "up" vector (and calculation intermediates) */ - //int done_debug = 0; + double nx, ny, nz; /* "normal" vector */ + double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda */ + double ux, uy, uz; /* "up" vector */ + double rx, ry, rz; /* "right" vector */ + int nrr = 0; + int n = 0; flist = image_feature_list_new(); @@ -67,15 +73,12 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * kx = nx / image->lambda; ky = ny / image->lambda; kz = nz / image->lambda; /* This is the centre of the Ewald sphere */ - //reflection_add(ctx->reflectionlist, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1); - /* Determine where "up" is - * See above. */ - uxt = 0.0; uyt = 1.0; uzt = 0.0; - ux = uxt; uy = cos(tilt)*uyt + sin(tilt)*uzt; uz = -sin(tilt)*uyt + cos(tilt)*uzt; - uxt = ux; uyt = uy; uzt = uz; - ux = uxt*cos(-omega) + uyt*-sin(omega); uy = -uxt*sin(omega) + uyt*cos(omega); uz = uzt; - //reflection_add(ctx->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2); + /* Determine where "up" is */ + mapping_rotate(0.0, 1.0, 0.0, &ux, &uy, &uz, omega, tilt); + + /* Determine where "right" is */ + mapping_rotate(1.0, 0.0, 0.0, &rx, &ry, &rz, omega, tilt); reflection = reflectionlist->reflections; do { @@ -92,7 +95,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * g_sq = modulus_squared(xl, yl, zl); gn = xl*nx + yl*ny + zl*nz; - /* Next, solve the relrod equation to calculate the excitation error */ + /* Next, solve the relrod equation to calculate + * the excitation error */ a = 1.0; b = 2.0*(gn - k); c = -2.0*gn*k + g_sq; @@ -106,10 +110,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * double xi, yi, zi; double gx, gy, gz; - double cx, cy, cz; double theta; double x, y; - double rx, ry, rz; /* Determine the intersection point */ xi = xl + s*nx; yi = yl + s*ny; zi = zl + s*nz; @@ -117,64 +119,36 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * /* Calculate Bragg angle */ gx = xi - kx; gy = yi - ky; - gz = zi - kz; /* This is the vector from the centre of the sphere to the intersection */ + gz = zi - kz; /* This is the vector from the centre of + * the sphere to the intersection */ theta = angle_between(-kx, -ky, -kz, gx, gy, gz); if ( theta > 0.0 ) { - double psi, disc; - - //reflection_add(ctx->reflectionlist, xl, yl, zl, 1, REFLECTION_GENERATED); - //reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_MARKER); - - /* Calculate azimuth of point in image (clockwise from "up", will be changed later) */ - cx = yi*nz-zi*ny; cy = nx*zi-nz*xi; cz = ny*xi-nx*yi; /* c = i x n */ - psi = angle_between(cx, cy, cz, ux, uy, uz); - - /* Finally, resolve the +/- Pi ambiguity from the previous step */ - rx = cy*nz-cz*ny; ry = nx*cz-nz*cx; rz = ny*cx-nx*cy; /* r = [i x n] x n */ - disc = angle_between(rx, ry, rz, ux, uy, uz); - // if ( (i==20) && !done_debug ) { - // reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_VECTOR_MARKER_3); - // reflection_add(ctx->reflectionlist, cx, cy, cz, 1, REFLECTION_VECTOR_MARKER_4); - // reflection_add(ctx->reflectionlist, rx, ry, rz, 1, REFLECTION_VECTOR_MARKER_4); - // printf("psi=%f deg, disc=%f deg\n", rad2deg(psi), rad2deg(disc)); - // } - if ( (psi >= M_PI_2) && (disc >= M_PI_2) ) { - psi -= M_PI_2; /* Case 1 */ - // if ( (i==20) && !done_debug ) printf("case 1\n"); - } else if ( (psi >= M_PI_2) && (disc < M_PI_2) ) { - psi = 3*M_PI_2 - psi; /* Case 2 */ - // if ( (i==20) && !done_debug ) printf("case 2\n"); - } else if ( (psi < M_PI_2) && (disc < M_PI_2) ) { - psi = 3*M_PI_2 - psi; /* Case 3 */ - // if ( (i==20) && !done_debug ) printf("case 3\n"); - } else if ( (psi < M_PI_2) && (disc >= M_PI_2) ) { - psi = 3*M_PI_2 + psi; /* Case 4 */ - // if ( (i==20) && !done_debug ) printf("case 4\n"); - } + double dx, dy, psi; - // if ( (i==20) && !done_debug ) printf("final psi=%f clockwise from 'up'\n", rad2deg(psi)); - psi = M_PI_2 - psi; /* Anticlockwise from "+x" instead of clockwise from "up" */ - // if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +x\n", rad2deg(psi)); + /* Calculate azimuth of point in image + * (anticlockwise from +x) */ + dx = xi*rx + yi*ry + zi*rz; + dy = xi*ux + yi*uy + zi*uz; + psi = atan2(dy, dx); - psi += omega; - // if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +tilt axis\n", rad2deg(psi)); - // if ( (i==20) && !done_debug ) done_debug = 1; - - /* Calculate image coordinates from polar representation */ + /* Get image coordinates from polar + * representation */ if ( image->fmode == FORMULATION_CLEN ) { x = image->camera_len*tan(theta)*cos(psi); y = image->camera_len*tan(theta)*sin(psi); x *= image->resolution; y *= image->resolution; - } else if ( image->fmode == FORMULATION_PIXELSIZE ) { + } else if ( image->fmode==FORMULATION_PIXELSIZE ) { x = tan(theta)*cos(psi) / image->lambda; y = tan(theta)*sin(psi) / image->lambda; x /= image->pixel_size; y /= image->pixel_size; } else { - fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n"); + fprintf(stderr, + "Unrecognised formulation mode " + "in reproject_get_reflections\n"); return NULL; } @@ -182,31 +156,45 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * y += image->y_centre; /* Sanity check */ - if ( (x>=0) && (xwidth) && (y>=0) && (yheight) ) { - - /* Record the reflection */ - image_add_feature_reflection(flist, x, y, image, reflection->intensity, reflection); - /* Intensity should be multiplied by relrod spike function except - * reprojected reflections aren't used quantitatively for anything */ - - //printf("Reflection at %f, %f\n", x, y); + if ( (x>=0) && (xwidth) + && (y>=0) && (yheight) ) { + + /* Record the reflection + * Intensity should be multiplied by + * relrod spike function except + * reprojected reflections aren't used + * quantitatively for anything + */ + image_add_feature_reflection(flist, x, + y, image, reflection->intensity, + reflection); + printf("adding %i %i %i\n", + reflection->h, + reflection->k, + reflection->l); + n++; } /* else it's outside the picture somewhere */ - } /* else this is the central beam so don't worry about it */ + } /* else this is the central beam */ } reflection = reflection->next; + nrr++; + } while ( reflection ); - /* Partner features only if the image has a feature list. This allows the test - * program to use this function to generate simulated data. */ + /* Partner features only if the image has a feature list. This allows + * the test programs to use this function to generate simulated data + */ if ( image->features != NULL ) { reproject_partner_features(flist, image); } + printf("processed %i, found %i\n", nrr, n); + return flist; } @@ -226,8 +214,20 @@ void reproject_cell_to_lattice(ControlContext *ctx) { /* Invalidate all the reprojected feature lists */ for ( i=0; iimages->n_images; i++ ) { - image_feature_list_free(ctx->images->images[i].rflist); - ctx->images->images[i].rflist = NULL; + + ImageRecord *image; + int j; + + image = &ctx->images->images[i]; + if ( image->rflist != NULL ) { + image_feature_list_free(image->rflist); + image->rflist = NULL; + } + + for ( j=0; jfeatures->n_features; j++ ) { + image->features->features[i].partner = NULL; + } + } } @@ -235,7 +235,65 @@ void reproject_cell_to_lattice(ControlContext *ctx) { /* Notify that ctx->cell has changed (also need to call displaywindow_update) */ void reproject_lattice_changed(ControlContext *ctx) { + int i; + int total_reprojected = 0; + int total_measured = 0; + int partnered_reprojected = 0; + int partnered_measured = 0; + reproject_cell_to_lattice(ctx); - displaywindow_update_imagestack(ctx->dw); + printf("%i reflections\n", ctx->cell_lattice->n_reflections); + for ( i=0; iimages->n_images; i++ ) { + + ImageRecord *image; + int j; + + image = &ctx->images->images[i]; + + /* Perform relrod projection */ + image->rflist = reproject_get_reflections(image, + ctx->cell_lattice); + + /* Loop over reprojected reflections */ + for ( j=0; jrflist->n_features; j++ ) { + + double d; + ImageFeature f; + ImageFeature *p; + + f = image->rflist->features[j]; + p = image->rflist->features[j].partner; + + if ( p != NULL ) { + d = distance(f.x, f.y, p->x, p->y); + partnered_reprojected++; + } + } + total_reprojected += image->rflist->n_features; + + /* Loop over measured reflections */ + for ( j=0; jfeatures->n_features; j++ ) { + + double d; + ImageFeature f; + ImageFeature *p; + + f = image->features->features[j]; + p = image->features->features[j].partner; + + if ( p != NULL ) { + d = distance(f.x, f.y, p->x, p->y); + partnered_measured++; + } + } + total_measured += image->features->n_features; + } + printf("%i images\n", ctx->images->n_images); + printf("%i measured reflections\n", total_measured); + printf("%i reprojected reflections\n", total_reprojected); + printf("%i partnered measured reflections\n", partnered_measured); + printf("%i partnered reprojected reflections\n", partnered_reprojected); + + displaywindow_update_imagestack(ctx->dw); } -- cgit v1.2.3