From 4262a193f2eea72dfd1a268a312f2ce5be71a6f1 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 31 Mar 2009 14:37:05 +0100 Subject: Reprojection improvements (previously on simple-search, lost, and now cherry-picked into master) New line-sphere intersection mathematics (Slightly) more numerically stable quadratic formula Conflicts: src/reproject.c --- src/reproject.c | 31 +++++++++++++++---------------- src/utils.c | 4 ++++ src/utils.h | 1 + 3 files changed, 20 insertions(+), 16 deletions(-) diff --git a/src/reproject.c b/src/reproject.c index c41aa3f..df85f4f 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -49,9 +49,9 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * ImageFeatureList *flist; Reflection *reflection; - double smax = 0.1e9; - double tilt, omega; - double nx, ny, nz, nxt, nyt, nzt; /* "normal" vector (and calculation intermediates) */ + double smax = 0.4e9; + 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; @@ -60,13 +60,10 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * 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 */ @@ -85,21 +82,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); -- cgit v1.2.3