Reprojection improvements

New line-sphere intersection mathematics
(Slightly) more numerically stable quadratic formula

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