diff options
author | Thomas White <taw27@cam.ac.uk> | 2008-11-04 17:41:25 +0000 |
---|---|---|
committer | Thomas White <taw27@cam.ac.uk> | 2008-11-04 17:41:25 +0000 |
commit | b44df6cd088c8cd1888c81d81f60502c9c5b588d (patch) | |
tree | 5489218a5727c99e3810a35dafc0273df2b00f3b | |
parent | cb76590c2cb5370aba709cee8a2e0514a9f423ff (diff) |
Reprojection improvements
New line-sphere intersection mathematics
(Slightly) more numerically stable quadratic formula
-rw-r--r-- | src/reproject.c | 19 | ||||
-rw-r--r-- | src/utils.c | 4 | ||||
-rw-r--r-- | src/utils.h | 1 |
3 files changed, 16 insertions, 8 deletions
diff --git a/src/reproject.c b/src/reproject.c index c594bab..48f60da 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -52,7 +52,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * ImageFeatureList *flist; Reflection *reflection; double smax = 0.1e9; - double tilt, omega; + double tilt, omega, k; double nx, ny, nz; /* "normal" vector */ double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda */ double ux, uy, uz; /* "up" vector */ @@ -62,6 +62,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * tilt = image->tilt; omega = image->omega; + k = 1.0/image->lambda; /* Calculate the (normalised) incident electron wavevector */ mapping_rotate(0.0, 0.0, 1.0, &nx, &ny, &nz, omega, tilt); @@ -80,21 +81,23 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * 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 */ diff --git a/src/utils.c b/src/utils.c index 30a6e34..01b14d9 100644 --- a/src/utils.c +++ b/src/utils.c @@ -39,6 +39,10 @@ double modulus(double x, double y, double z) { return sqrt(x*x + y*y + z*z); } +double modulus_squared(double x, double y, double z) { + return x*x + y*y + z*z; +} + double distance3d(double x1, double y1, double z1, double x2, double y2, double z2) { return modulus(x1-x2, y1-y2, z1-z2); } diff --git a/src/utils.h b/src/utils.h index 5c2e586..29ebf01 100644 --- a/src/utils.h +++ b/src/utils.h @@ -23,6 +23,7 @@ extern unsigned int biggest(signed int a, signed int b); extern unsigned int smallest(signed int a, signed int b); extern double distance(double x1, double y1, double x2, double y2); extern double modulus(double x, double y, double z); +extern double modulus_squared(double x, double y, double z); extern double angle_between(double x1, double y1, double z1, double x2, double y2, double z2); extern double angle_between_d(double x1, double y1, double z1, double x2, double y2, double z2); extern double lambda(double voltage); |