From b2a198a4a7935d4c81c0b7044d9f89e3c6932472 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 20 May 2014 17:43:18 +0200 Subject: Add Gaussian partiality model --- libcrystfel/src/geometry.c | 42 +++++++++++++++++++++++--------- libcrystfel/src/geometry.h | 8 ++++--- src/hrs-scaling.c | 57 +++----------------------------------------- src/partialator.c | 2 ++ src/post-refinement.c | 56 +++++++++++++++++++++++++++---------------- tests/pr_l_gradient_check.c | 12 ++++++++-- tests/pr_p_gradient_check.c | 12 ++++++++-- tests/pr_pl_gradient_check.c | 12 ++++++++-- 8 files changed, 107 insertions(+), 94 deletions(-) diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 0c0a09b8..357f64f8 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -3,11 +3,11 @@ * * Geometry of diffraction * - * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2013 Thomas White + * 2010-2014 Thomas White * * This file is part of CrystFEL. * @@ -34,6 +34,7 @@ #include #include #include +#include #include "utils.h" #include "cell.h" @@ -102,24 +103,44 @@ static signed int locate_peak(double x, double y, double z, double k, } -static double partiality(double rlow, double rhigh, double r) +static double partiality(PartialityModel pmodel, double rlow, double rhigh, + double r) { double qlow, qhigh; double plow, phigh; + const double ng = 2.6; /* Calculate degrees of penetration */ qlow = (rlow + r)/(2.0*r); qhigh = (rhigh + r)/(2.0*r); /* Convert to partiality */ - plow = 3.0*pow(qlow,2.0) - 2.0*pow(qlow,3.0); - phigh = 3.0*pow(qhigh,2.0) - 2.0*pow(qhigh,3.0); + switch ( pmodel ) { + + default: + case PMODEL_UNITY: + plow = 1.0; + phigh = 0.0; + break; + + case PMODEL_SPHERE: + plow = 3.0*pow(qlow,2.0) - 2.0*pow(qlow,3.0); + phigh = 3.0*pow(qhigh,2.0) - 2.0*pow(qhigh,3.0); + break; + + case PMODEL_GAUSSIAN: + plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*r)); + phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*r)); + break; + + } return plow - phigh; } static Reflection *check_reflection(struct image *image, Crystal *cryst, + PartialityModel pmodel, signed int h, signed int k, signed int l, double xl, double yl, double zl) { @@ -204,7 +225,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, } /* Calculate partiality */ - part = partiality(rlow, rhigh, pr); + part = partiality(pmodel, rlow, rhigh, pr); /* Add peak to list */ refl = reflection_new(h, k, l); @@ -298,7 +319,8 @@ RefList *find_intersections(struct image *image, Crystal *cryst) yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - refl = check_reflection(image, cryst, h, k, l, xl, yl, zl); + refl = check_reflection(image, cryst, PMODEL_SPHERE, + h, k, l, xl, yl, zl); if ( refl != NULL ) { add_refl_to_list(refl, reflections); @@ -401,9 +423,6 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel, int n = 0; if ( pmodel == PMODEL_UNITY ) { - /* It isn't strictly necessary to set the partialities to 1, - * because the scaling stuff will just a correction factor of - * 1 anyway. This is just to help things make sense. */ set_unity_partialities(cryst); return; } @@ -430,7 +449,8 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel, yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - vals = check_reflection(image, cryst, h, k, l, xl, yl, zl); + vals = check_reflection(image, cryst, pmodel, + h, k, l, xl, yl, zl); if ( vals == NULL ) { diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index fa7ed53b..7188f6ed 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -3,12 +3,12 @@ * * Geometry of diffraction * - * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * * Authors: - * 2010-2013 Thomas White + * 2010-2014 Thomas White * 2012 Richard Kirian * * This file is part of CrystFEL. @@ -49,6 +49,7 @@ extern "C" { * @PMODEL_SPHERE : Intersection of sphere with excited volume of reciprocal * space. * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1. + * @PMODEL_GAUSSIAN : Gaussian profiles in 3D * * A %PartialityModel describes a geometrical model which can be used to * calculate spot partialities and Lorentz correction factors. @@ -57,6 +58,7 @@ typedef enum { PMODEL_SPHERE, PMODEL_UNITY, + PMODEL_GAUSSIAN, } PartialityModel; diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 76d1ff75..ac9091ed 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -119,24 +119,7 @@ static void run_scale_job(void *vwargs, int cookie) Ih = get_intensity(r); } - /* If you change this, be sure to also change - * run_merge_job() and run_esd_job(). */ - switch ( wargs->pmodel ) { - - case PMODEL_UNITY : - corr = 1.0; - break; - - case PMODEL_SPHERE : - corr = get_partiality(refl) * get_lorentz(refl); - break; - - default : - ERROR("Unrecognised partiality model!\n"); - abort(); - break; - - } + corr = get_partiality(refl) * get_lorentz(refl); Ihl = get_intensity(refl) / corr; @@ -285,24 +268,7 @@ static void run_merge_job(void *vwargs, int cookie) } - /* If you change this, be sure to also change - * run_scale_job() and run_esd_job(). */ - switch ( wargs->pmodel ) { - - case PMODEL_UNITY : - corr = 1.0; - break; - - case PMODEL_SPHERE : - corr = get_partiality(refl) * get_lorentz(refl); - break; - - default : - ERROR("Unrecognised partiality model!\n"); - abort(); - break; - - } + corr = get_partiality(refl) * get_lorentz(refl); Ihl = get_intensity(refl) / corr; @@ -421,24 +387,7 @@ static void run_esd_job(void *vwargs, int cookie) num = get_temp1(f); - /* If you change this, be sure to also change - * run_scale_job() and run_merge_job(). */ - switch ( wargs->pmodel ) { - - case PMODEL_UNITY : - corr = 1.0; - break; - - case PMODEL_SPHERE : - corr = get_partiality(refl) * get_lorentz(refl); - break;; - - default : - ERROR("Unrecognised partiality model!\n"); - abort(); - break; - - } + corr = get_partiality(refl) * get_lorentz(refl); Ih = get_intensity(f); Ihl = G * get_intensity(refl) / corr; diff --git a/src/partialator.c b/src/partialator.c index bc8704b0..c1dedad8 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -468,6 +468,8 @@ int main(int argc, char *argv[]) pmodel = PMODEL_SPHERE; } else if ( strcmp(pmodel_str, "unity") == 0 ) { pmodel = PMODEL_UNITY; + } else if ( strcmp(pmodel_str, "gaussian") == 0 ) { + pmodel = PMODEL_GAUSSIAN; } else { ERROR("Unknown partiality model '%s'.\n", pmodel_str); return 1; diff --git a/src/post-refinement.c b/src/post-refinement.c index d5d8a4d1..8fd29783 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2013 Thomas White + * 2010-2014 Thomas White * * This file is part of CrystFEL. * @@ -52,39 +52,55 @@ #define MAX_CYCLES (10) -/* Returns dp/dr at "r" */ -static double partiality_gradient(double r, double profile_radius) +static double dpdq(double r, double profile_radius, PartialityModel pmodel) { - double q, dpdq, dqdr; + double q; + double ng = 3.0; /* Calculate degree of penetration */ q = (r + profile_radius)/(2.0*profile_radius); /* dp/dq */ - dpdq = 6.0*(q-pow(q, 2.0)); + switch ( pmodel ) { + + default: + case PMODEL_UNITY: + return 0.0; + + case PMODEL_SPHERE: + return 6.0*(q-pow(q, 2.0)); + + case PMODEL_GAUSSIAN: + /* The flat sphere model is close enough */ + return 6.0*(q-pow(q, 2.0)); + + } +} + + +/* Returns dp/dr at "r" */ +static double partiality_gradient(double r, double profile_radius, + PartialityModel pmodel) +{ + double dqdr; /* dq/dr */ dqdr = 1.0 / (2.0*profile_radius); - return dpdq * dqdr; + return dpdq(r, profile_radius, pmodel) * dqdr; } /* Returns dp/drad at "r" */ -static double partiality_rgradient(double r, double profile_radius) +static double partiality_rgradient(double r, double profile_radius, + PartialityModel pmodel) { - double q, dpdq, dqdrad; - - /* Calculate degree of penetration */ - q = (r + profile_radius)/(2.0*profile_radius); - - /* dp/dq */ - dpdq = 6.0*(q-pow(q, 2.0)); + double dqdrad; /* dq/drad */ dqdrad = -0.5 * r * pow(profile_radius, -2.0); - return dpdq * dqdrad; + return dpdq(r, profile_radius, pmodel) * dqdrad; } @@ -144,12 +160,12 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) /* Calculate the gradient of partiality wrt excitation error. */ if ( clamp_low == 0 ) { - glow = partiality_gradient(rlow, r); + glow = partiality_gradient(rlow, r, pmodel); } else { glow = 0.0; } if ( clamp_high == 0 ) { - ghigh = partiality_gradient(rhigh, r); + ghigh = partiality_gradient(rhigh, r, pmodel); } else { ghigh = 0.0; } @@ -163,8 +179,8 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) return (ds*glow + ds*ghigh) / 2.0; case REF_R : - gr = partiality_rgradient(rlow, r); - gr -= partiality_rgradient(rhigh, r); + gr = partiality_rgradient(rlow, r, pmodel); + gr -= partiality_rgradient(rhigh, r, pmodel); return gr; /* Cell parameters and orientation */ diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c index 4b3894ba..72d353f9 100644 --- a/tests/pr_l_gradient_check.c +++ b/tests/pr_l_gradient_check.c @@ -279,7 +279,6 @@ int main(int argc, char *argv[]) Crystal *cr; struct quaternion orientation; int i; - const PartialityModel pmodel = PMODEL_SPHERE; int fail = 0; int quiet = 0; int plot = 0; @@ -337,10 +336,19 @@ int main(int argc, char *argv[]) rng = gsl_rng_alloc(gsl_rng_mt19937); - for ( i=0; i<1; i++ ) { + for ( i=0; i<2; i++ ) { UnitCell *rot; double val; + PartialityModel pmodel; + + if ( i == 0 ) { + pmodel = PMODEL_SPHERE; + STATUS("Testing flat sphere model:\n"); + } else { + pmodel = PMODEL_GAUSSIAN; + STATUS("Testing Gaussian model:\n"); + } orientation = random_quaternion(rng); rot = cell_rotate(cell, orientation); diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c index c0582036..c3b55cfa 100644 --- a/tests/pr_p_gradient_check.c +++ b/tests/pr_p_gradient_check.c @@ -392,7 +392,6 @@ int main(int argc, char *argv[]) Crystal *cr; struct quaternion orientation; int i; - const PartialityModel pmodel = PMODEL_SPHERE; int fail = 0; int quiet = 0; int plot = 0; @@ -450,10 +449,19 @@ int main(int argc, char *argv[]) rng = gsl_rng_alloc(gsl_rng_mt19937); - for ( i=0; i<1; i++ ) { + for ( i=0; i<2; i++ ) { UnitCell *rot; double val; + PartialityModel pmodel; + + if ( i == 0 ) { + pmodel = PMODEL_SPHERE; + STATUS("Testing flat sphere model:\n"); + } else { + pmodel = PMODEL_GAUSSIAN; + STATUS("Testing Gaussian model:\n"); + } orientation = random_quaternion(rng); rot = cell_rotate(cell, orientation); diff --git a/tests/pr_pl_gradient_check.c b/tests/pr_pl_gradient_check.c index 785d9973..ddbc1840 100644 --- a/tests/pr_pl_gradient_check.c +++ b/tests/pr_pl_gradient_check.c @@ -394,7 +394,6 @@ int main(int argc, char *argv[]) Crystal *cr; struct quaternion orientation; int i; - const PartialityModel pmodel = PMODEL_SPHERE; int fail = 0; int quiet = 0; int plot = 0; @@ -452,10 +451,19 @@ int main(int argc, char *argv[]) rng = gsl_rng_alloc(gsl_rng_mt19937); - for ( i=0; i<1; i++ ) { + for ( i=0; i<2; i++ ) { UnitCell *rot; double val; + PartialityModel pmodel; + + if ( i == 0 ) { + pmodel = PMODEL_SPHERE; + STATUS("Testing flat sphere model:\n"); + } else { + pmodel = PMODEL_GAUSSIAN; + STATUS("Testing Gaussian model:\n"); + } orientation = random_quaternion(rng); rot = cell_rotate(cell, orientation); -- cgit v1.2.3