From 363d209ab72bd42d401d86ce9137f9392f4b8e19 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 15 Oct 2008 18:35:26 +0100 Subject: Use dot product projection for reprojection The refinement code uses a projection based on taking dot products with the x and y directions. Do the same thing for reprojection to find the azimuth of the scattered beam, instead of the NASTY stuff currently used. --- src/reproject.c | 32 ++++++++++---------------------- 1 file changed, 10 insertions(+), 22 deletions(-) diff --git a/src/reproject.c b/src/reproject.c index 46d59c9..e145186 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -55,7 +55,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * 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 */ + flist = image_feature_list_new(); tilt = image->tilt; @@ -72,6 +73,9 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * mapping_rotate(0.0, 1.0, 0.0, &ux, &uy, &uz, omega, tilt); reflection_add(ctx->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2); + /* Determine where "right" is */ + mapping_rotate(1.0, 0.0, 0.0, &rx, &ry, &rz, omega, tilt); + reflection = reflectionlist->reflections; do { @@ -99,10 +103,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; @@ -115,26 +117,12 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * if ( theta > 0.0 ) { - double psi, disc; - - /* 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 ( (psi >= M_PI_2) && (disc >= M_PI_2) ) { - psi -= M_PI_2; /* Case 1 */ - } else if ( (psi >= M_PI_2) && (disc < M_PI_2) ) { - psi = 3*M_PI_2 - psi; /* Case 2 */ - } else if ( (psi < M_PI_2) && (disc < M_PI_2) ) { - psi = 3*M_PI_2 - psi; /* Case 3 */ - } else if ( (psi < M_PI_2) && (disc >= M_PI_2) ) { - psi = 3*M_PI_2 + psi; /* Case 4 */ - } + double dx, dy, psi; - psi = M_PI_2 - psi; /* Anticlockwise from "+x" instead of clockwise from "up" */ + /* 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); /* Calculate image coordinates from polar representation */ if ( image->fmode == FORMULATION_CLEN ) { -- cgit v1.2.3