diff options
| author | Thomas White <taw27@cam.ac.uk> | 2008-10-15 18:35:26 +0100 |
|---|---|---|
| committer | Thomas White <taw27@cam.ac.uk> | 2008-10-15 18:35:26 +0100 |
| commit | 363d209ab72bd42d401d86ce9137f9392f4b8e19 (patch) | |
| tree | d94b7ddb6abedebfc54cb883a2917a27e155503d | |
| parent | a066052be7beab1a2657dced4a03e8be4d9e3bd6 (diff) | |
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.
| -rw-r--r-- | src/reproject.c | 32 |
1 files 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 ) { |