aboutsummaryrefslogtreecommitdiff
path: root/src/reproject.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/reproject.c')
-rw-r--r--src/reproject.c113
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 */
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-01 16:30:53 +0000