diff options
Diffstat (limited to 'src/reproject.c')
-rw-r--r-- | src/reproject.c | 113 |
1 files changed, 38 insertions, 75 deletions
diff --git a/src/reproject.c b/src/reproject.c index e1ce13e..48f60da 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) { @@ -32,7 +33,8 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) { ImageFeature *partner; int idx; - partner = image_feature_closest(image->features, rflist->features[i].x, rflist->features[i].y, &d, &idx); + partner = image_feature_closest(image->features, rflist->features[i].x, rflist->features[i].y, + &d, &idx); if ( (d <= 20.0) && partner ) { rflist->features[i].partner = partner; @@ -45,61 +47,57 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) { } -ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist) { +ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist, ControlContext *ctx) { ImageFeatureList *flist; Reflection *reflection; double smax = 0.1e9; - double tilt, omega; - double nx, ny, nz, nxt, nyt, nzt; /* "normal" vector (and calculation intermediates) */ + 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 ux, uy, uz; /* "up" vector */ + double rx, ry, rz; /* "right" vector */ + flist = image_feature_list_new(); tilt = image->tilt; omega = image->omega; + k = 1.0/image->lambda; - /* Calculate the (normalised) incident electron wavevector - * n is rotated "into" the reconstruction, so only one omega step. */ - nxt = 0.0; nyt = 0.0; nzt = 1.0; - nx = nxt; ny = cos(tilt)*nyt + sin(tilt)*nzt; nz = -sin(tilt)*nyt + cos(tilt)*nzt; - nxt = nx; nyt = ny; nzt = nz; - nx = nxt*cos(-omega) + nyt*sin(-omega); ny = -nxt*sin(-omega) + nyt*cos(-omega); nz = nzt; + /* Calculate the (normalised) incident electron wavevector */ + mapping_rotate(0.0, 0.0, 1.0, &nx, &ny, &nz, omega, tilt); 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 { double xl, yl, zl; double a, b, c; - double A1, A2, s1, s2, s; + double s1, s2, s, t; + double g_sq, gn; /* Get the coordinates of the reciprocal lattice point */ xl = reflection->x; yl = reflection->y; zl = reflection->z; + g_sq = modulus_squared(xl, yl, zl); + gn = xl*nx + yl*ny + zl*nz; /* Next, solve the relrod equation to calculate the excitation error */ a = 1.0; - b = -2.0*(xl*nx + yl*ny + zl*nz); - c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda); - A1 = (-b + sqrt(b*b-4.0*a*c))/(2.0*a); - A2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a); - s1 = 1.0/image->lambda - A1; - s2 = 1.0/image->lambda - A2; + b = 2.0*(gn - k); + c = -2.0*gn*k + g_sq; + t = -0.5*(b + sign(b)*sqrt(b*b - 4.0*a*c)); + s1 = t/a; + s2 = c/t; if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2; /* Skip this reflection if s is large */ @@ -107,10 +105,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; @@ -123,45 +119,12 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * if ( theta > 0.0 ) { - double psi, disc; + double dx, dy, psi; - //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"); - } - - // 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)); - - 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 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 ) { @@ -175,7 +138,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * 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; } @@ -185,13 +148,13 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * /* Sanity check */ if ( (x>=0) && (x<image->width) && (y>=0) && (y<image->height) ) { - /* 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 */ + /* 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("Reflection at %f, %f\n", x, y); - } /* else it's outside the picture somewhere */ } /* else this is the central beam so don't worry about it */ |