diff options
Diffstat (limited to 'libcrystfel')
| -rw-r--r-- | libcrystfel/src/geometry.c | 133 | ||||
| -rw-r--r-- | libcrystfel/src/peaks.c | 24 | ||||
| -rw-r--r-- | libcrystfel/src/utils.h | 23 |
3 files changed, 151 insertions, 29 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index e5f33d64..4d25b89f 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -239,39 +239,135 @@ static double random_partiality(signed int h, signed int k, signed int l, } +static inline double safe_khalf(const double xl, + const double yl, + const double zl) +{ + if ( zl > 0.0 ) return 0.0/0.0; + return -(xl*xl+yl*yl+zl*zl) / (2.0*zl); +} + + static Reflection *check_reflection(struct image *image, Crystal *cryst, signed int h, signed int k, signed int l, double xl, double yl, double zl, Reflection *updateme) { Reflection *refl; - double R, top; - double kmin, kmax, k0, knom, k1, khalf; + double R; + int i, n; + double knom, khalf; double dcs, exerr; + double partiality, mean_kpred, M2_kpred; + double sumw_k, mean_k, M2_k; + + /* This arbitrary value is there to mimic previous behaviour */ + const double min_partiality = exp(-0.5*1.7*1.7); /* Don't predict 000 */ if ( (updateme == NULL) && (abs(h)+abs(k)+abs(l) == 0) ) return NULL; - /* Calculate the limiting wavelengths, lambda0 and lambda1 - * = 1/k0 and 1/k1 respectively */ R = fabs(crystal_get_profile_radius(cryst)); - top = R*R - xl*xl - yl*yl - zl*zl; - k0 = top/(2.0*(zl+R)); - khalf = (- xl*xl - yl*yl - zl*zl) / (2.0*zl); - k1 = top/(2.0*(zl-R)); + n = spectrum_get_num_gaussians(image->spectrum); + assert(n > 0); + + partiality = 0.0; + mean_kpred = 0.0; + M2_kpred = 0.0; + sumw_k = 0.0; + mean_k = 0.0; + M2_k = 0.0; + for ( i=0; i<n; i++ ) { + + struct gaussian g; + double kpred; + double exerr2, x, y, z, norm; + double sigma_proj, w0, w1; + double sigma2, exponent, overlap_integral; + + g = spectrum_get_gaussian(image->spectrum, i); + + /* Project lattice point onto Ewald sphere */ + x = xl; + y = yl; + z = zl + g.kcen; + norm = 1.0/sqrt(x*x+y*y+z*z); + x *= norm; + y *= norm; + z *= norm; + z -= 1.0; + + /* Width of Ewald sphere in the direction of the projection */ + sigma_proj = sqrt(x*x+y*y+z*z)*g.sigma; + + mean_variance(g.kcen, g.area, &sumw_k, &mean_k, &M2_k); + M2_k += g.area * g.sigma * g.sigma; + w0 = 1.0/(R*R); + w1 = 1.0/(sigma_proj*sigma_proj); + + x *= g.kcen; + y *= g.kcen; + z *= g.kcen; + + /* Three because the general case fails in extreme cases */ + if ( w0 / w1 <= DBL_MIN ) { + + /* 'Laue' corner case */ + kpred = g.kcen; + exerr2 = g.kcen - safe_khalf(xl, yl, zl); + exerr2 *= exerr2; + + } else if ( w1 / w0 <= DBL_MIN ) { + + /* 'Monochromatic' corner case */ + kpred = safe_khalf(xl,yl,zl); + exerr2 = g.kcen - kpred; + exerr2*= exerr2; - /* The reflection is excited if any of the reflection is within 2sigma - * of the nominal * wavelength of the X-ray beam - * (NB image->bw is full width) */ - kmin = 1.0/(image->lambda + image->lambda*image->bw); - knom = 1.0/image->lambda; - kmax = 1.0/(image->lambda - image->lambda*image->bw); - if ( (updateme == NULL) && ((k1>kmax) || (k0<kmin)) ) return NULL; + } else { + + /* General case */ + + /* Closest point on Ewald sphere. + * Project zl to 0, bit of a hack... */ + const double zlp0 = zl>0?zl:0; + exerr2 = (x-xl)*(x-xl) + (y-yl)*(y-yl) + (z-zl)*(z-zl); + + /* Weighted average between projected lattice point + * and Ewald sphere */ + x = ( xl *w0 + x*w1 ) / ( w0 + w1 ); + y = ( yl *w0 + y*w1 ) / ( w0 + w1 ); + z = ( zlp0*w0 + z*w1 ) / ( w0 + w1 ); + kpred = safe_khalf(x,y,z); + + } + sigma2 = R*R + sigma_proj*sigma_proj; + exponent = - 0.5 * exerr2 / sigma2; + if ( exponent > -700.0 ) { + overlap_integral = exp(exponent) * sqrt(2*M_PI*R*R) / sqrt(2*M_PI*sigma2); + } else { + exponent = 0.0; + } + + mean_variance(kpred, g.area*overlap_integral, + &partiality, &mean_kpred, &M2_kpred); + + } + + /* Revert the 'Lorentz' factor */ + partiality *= sqrt( ( R*R + M2_k/sumw_k) / ( R*R ) ); + + if ( (updateme == NULL) && ( partiality < min_partiality ) ) return NULL; /* Calculate excitation error */ + knom = 1.0/image->lambda; dcs = distance3d(0.0, 0.0, -knom, xl, yl, zl); exerr = 1.0/image->lambda - dcs; /* Loss of precision */ + /* Could also estimate excitation error like this, but it loses the sign + * (which is needed elsewhere): + * exerr = R * sqrt(-2*log(partiality)); */ + if ( updateme == NULL ) { refl = reflection_new(h, k, l); } else { @@ -284,7 +380,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, if ( (image->det != NULL) && (updateme != NULL) ) { double fs, ss; - locate_peak_on_panel(xl, yl, zl, knom, + locate_peak_on_panel(xl, yl, zl, mean_kpred, get_panel(updateme), &fs, &ss); set_detector_pos(refl, fs, ss); @@ -296,7 +392,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, double fs, ss; /* Position on detector */ signed int p; /* Panel number */ - p = locate_peak(xl, yl, zl, knom, + p = locate_peak(xl, yl, zl, mean_kpred, image->det, &fs, &ss); if ( p == -1 ) { reflection_free(refl); @@ -307,7 +403,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, } - set_kpred(refl, knom); + khalf = (- xl*xl - yl*yl - zl*zl) / (2.0*zl); + set_kpred(refl, mean_kpred); set_khalf(refl, khalf); set_exerr(refl, exerr); set_lorentz(refl, 1.0); diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 6da79e4f..f7433383 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -812,34 +812,36 @@ static int compare_double(const void *av, const void *bv) double estimate_peak_resolution(ImageFeatureList *peaks, double lambda) { - int i, n, j; + int i, npk, ncut; double *rns; double max_res; - n = image_feature_count(peaks); - rns = malloc(n*sizeof(double)); + npk = image_feature_count(peaks); + + /* No peaks -> no resolution! */ + if ( npk == 0 ) return 0.0; + + rns = malloc(npk*sizeof(double)); if ( rns == NULL ) return -1.0; /* Get resolution values for all peaks */ - j = 0; - for ( i=0; i<n; i++ ) { + for ( i=0; i<npk; i++ ) { struct imagefeature *f; struct rvec r; f = image_get_feature(peaks, i); - if ( f == NULL ) continue; r = get_q_for_panel(f->p, f->fs, f->ss, NULL, 1.0/lambda); - rns[j++] = modulus(r.u, r.v, r.w); + rns[i] = modulus(r.u, r.v, r.w); } /* Slightly horrible outlier removal */ - qsort(rns, j, sizeof(double), compare_double); - n = j/50; - if ( n < 2 ) n = 2; - max_res = rns[(j-1)-n]; + qsort(rns, npk, sizeof(double), compare_double); + ncut = npk/50; + if ( ncut < 2 ) ncut = 0; + max_res = rns[(npk-1)-ncut]; free(rns); return max_res; diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index cd21bd93..02aca4c7 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -35,6 +35,7 @@ #include <math.h> #include <complex.h> +#include <float.h> #include <string.h> #include <stdlib.h> #include <pthread.h> @@ -269,6 +270,28 @@ static inline struct quaternion invalid_quaternion(void) } +/** + * \param x value + * \param w weight + * \param sumw pointer to accumulator variable for the sum of weights + * \param mean pointer to online mean value + * \param M2 pointer to online variance times sumw + * + * Function to compute mean and variance stably + */ +static inline void mean_variance(const double x, const double w, + double *sumw, double *mean, double *M2) +{ + const double temp = w + *sumw; + const double delta = x - *mean; + const double R = delta * w / temp; + + if ( w < DBL_MIN ) return; + + *mean += R; + *M2 += *sumw*delta*R; + *sumw = temp; +} #ifdef __cplusplus } |