From 4ea397e11edfb5ef711e52a0d4088943df052d7b Mon Sep 17 00:00:00 2001 From: taw27 Date: Tue, 12 Aug 2008 16:30:06 +0000 Subject: Improve numerical stability of reprojection git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@285 bf6ca9ba-c028-0410-8290-897cf20841d1 --- src/reproject.c | 11 +++++++---- src/utils.c | 5 +++++ src/utils.h | 1 + 3 files changed, 13 insertions(+), 4 deletions(-) diff --git a/src/reproject.c b/src/reproject.c index e1ce13e..d071f9c 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -85,7 +85,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * double xl, yl, zl; double a, b, c; - double A1, A2, s1, s2, s; + double A1, A2, s1, s2, s, temp; /* Get the coordinates of the reciprocal lattice point */ xl = reflection->x; @@ -95,9 +95,12 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * /* 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); + c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda); /* FIXME: Don't think this is stable */ + printf("%e %e\n", xl*xl + yl*yl + zl*zl, 1.0/(image->lambda*image->lambda)); + /* Solve the quadratic equation */ + temp = -0.5 * (b + sign(b)*sqrt(b*b - 4*a*c)); + A1 = temp / a; + A2 = c / temp; s1 = 1.0/image->lambda - A1; s2 = 1.0/image->lambda - A2; if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2; diff --git a/src/utils.c b/src/utils.c index 3fe2e23..30a6e34 100644 --- a/src/utils.c +++ b/src/utils.c @@ -112,4 +112,9 @@ void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v, const char *pr } +int sign(double a) { + if ( a < 0 ) return -1; + if ( a > 0 ) return +1; + return 0; +} diff --git a/src/utils.h b/src/utils.h index 5736739..5c2e586 100644 --- a/src/utils.h +++ b/src/utils.h @@ -30,6 +30,7 @@ extern double distance3d(double x1, double y1, double z1, double x2, double y2, extern size_t skipspace(const char *s); extern void chomp(char *s); extern void matrix_vector_show(const gsl_matrix *m, const gsl_vector *v, const char *prefix); +extern int sign(double a); #define rad2deg(a) ((a)*180/M_PI) #define deg2rad(a) ((a)*M_PI/180) -- cgit v1.2.3