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

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);