aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/src/integration.c223
-rw-r--r--libcrystfel/src/integration.h6
-rw-r--r--libcrystfel/src/stream.c8
3 files changed, 121 insertions, 116 deletions
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index e753dea9..76f98802 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -177,7 +177,6 @@ struct intcontext
int n_saturated;
int n_implausible;
- double limit;
IntDiag int_diag;
signed int int_diag_h;
@@ -1530,7 +1529,6 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
struct intcontext *ic, UnitCell *cell)
{
double pfs, pss;
- signed int h, k, l;
struct peak_box *bx;
int pn;
struct panel *p;
@@ -1540,7 +1538,6 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
double intensity;
double sigma;
int saturated;
- double one_over_d;
int r;
double bgmean, sig2_bg, sig2_poisson, aduph;
@@ -1621,10 +1618,6 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
set_mean_bg(refl, bgmean);
set_peak(refl, bx->peak);
- get_indices(refl, &h, &k, &l);
- one_over_d = resolution(cell, h, k, l);
- if ( one_over_d > ic->limit ) ic->limit = one_over_d;
-
/* Update position */
pfs += bx->offs_fs;
pss += bx->offs_ss;
@@ -1639,6 +1632,94 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
}
+static int compare_double(const void *av, const void *bv)
+{
+ double a = *(double *)av;
+ double b = *(double *)bv;
+ if ( a > b ) return 1;
+ if ( a < b ) return -1;
+ return 0;
+}
+
+
+static double estimate_resolution(UnitCell *cell, ImageFeatureList *flist)
+{
+ int i;
+ const double min_dist = 0.25;
+ double max_res = 0.0;
+ double *acc;
+ int n_acc = 0;
+ int max_acc = 1024;
+ double av;
+ int n;
+
+ acc = malloc(max_acc*sizeof(double));
+ if ( acc == NULL ) {
+ ERROR("Allocation failed during estimate_resolution!\n");
+ return INFINITY;
+ }
+
+ for ( i=0; i<image_feature_count(flist); i++ ) {
+
+ struct imagefeature *f;
+ double h, k, l, hd, kd, ld;
+
+ /* Assume all image "features" are genuine peaks */
+ f = image_get_feature(flist, i);
+ if ( f == NULL ) continue;
+
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+
+ cell_get_cartesian(cell,
+ &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+ /* Decimal and fractional Miller indices of nearest
+ * reciprocal lattice point */
+ hd = f->rx * ax + f->ry * ay + f->rz * az;
+ kd = f->rx * bx + f->ry * by + f->rz * bz;
+ ld = f->rx * cx + f->ry * cy + f->rz * cz;
+ h = lrint(hd);
+ k = lrint(kd);
+ l = lrint(ld);
+
+ /* Check distance */
+ if ( (fabs(h - hd) < min_dist)
+ && (fabs(k - kd) < min_dist)
+ && (fabs(l - ld) < min_dist) )
+ {
+ double res = 2.0*resolution(cell, h, k, l); /* 1/d */
+ acc[n_acc++] = res;
+ if ( n_acc == max_acc ) {
+ max_acc += 1024;
+ acc = realloc(acc, max_acc*sizeof(double));
+ if ( acc == NULL ) {
+ ERROR("Allocation failed during"
+ " estimate_resolution!\n");
+ return INFINITY;
+ }
+ }
+ }
+
+ }
+
+ if ( n_acc < 3 ) {
+ STATUS("WARNING: Too few peaks to estimate resolution.\n");
+ return 0.0;
+ }
+
+ /* Slightly horrible outlier removal */
+ qsort(acc, n_acc, sizeof(double), compare_double);
+ n = n_acc/50;
+ if ( n < 2 ) n = 2;
+ max_res = acc[(n_acc-1)-n];
+
+ free(acc);
+ return max_res;
+}
+
+
static void integrate_rings(IntegrationMethod meth, Crystal *cr,
struct image *image, IntDiag int_diag,
signed int idh, signed int idk, signed int idl,
@@ -1670,7 +1751,6 @@ static void integrate_rings(IntegrationMethod meth, Crystal *cr,
ic.int_diag_h = idh;
ic.int_diag_k = idk;
ic.int_diag_l = idl;
- ic.limit = 0.0;
ic.meth = meth;
if ( init_intcontext(&ic) ) {
ERROR("Failed to initialise integration.\n");
@@ -1690,7 +1770,6 @@ static void integrate_rings(IntegrationMethod meth, Crystal *cr,
free_intcontext(&ic);
crystal_set_num_saturated_reflections(cr, ic.n_saturated);
- crystal_set_resolution_limit(cr, ic.limit);
crystal_set_reflections(cr, list);
if ( ic.n_implausible ) {
@@ -1700,113 +1779,24 @@ static void integrate_rings(IntegrationMethod meth, Crystal *cr,
}
-static void UNUSED resolution_cutoff(RefList *list, UnitCell *cell)
+static void apply_resolution_cutoff(Crystal *cr, double res)
{
- double *rmins;
- double *rmaxs;
- double rmin, rmax;
- double total_vol, vol_per_shell;
- int i;
- int *sh_nref;
- double *sh_snr;
- Reflection **highest_snr;
Reflection *refl;
RefListIterator *iter;
+ UnitCell *cell;
- const int nshells = 100;
- const double min_snr = 10.0;
-
- rmins = malloc(nshells*sizeof(double));
- rmaxs = malloc(nshells*sizeof(double));
- sh_nref = malloc(nshells*sizeof(int));
- sh_snr = malloc(nshells*sizeof(double));
- highest_snr = malloc(nshells*sizeof(Reflection *));
-
- if ( (rmins==NULL) || (rmaxs==NULL) || (sh_nref==NULL)
- || (sh_snr==NULL) || (highest_snr==NULL) )
- {
- ERROR("Couldn't allocate memory for resolution shells.\n");
- return;
- }
-
- resolution_limits(list, cell, &rmin, &rmax);
-
- total_vol = pow(rmax, 3.0) - pow(rmin, 3.0);
- vol_per_shell = total_vol / nshells;
- rmins[0] = rmin;
- for ( i=1; i<nshells; i++ ) {
-
- double r;
-
- r = vol_per_shell + pow(rmins[i-1], 3.0);
- r = pow(r, 1.0/3.0);
-
- /* Shells of constant volume */
- rmaxs[i-1] = r;
- rmins[i] = r;
-
- }
- rmaxs[nshells-1] = rmax;
-
- for ( i=0; i<nshells; i++ ) {
- sh_nref[i] = 0;
- sh_snr[i] = 0.0;
- }
+ cell = crystal_get_cell(cr);
- for ( refl = first_refl(list, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
signed int h, k, l;
- double res, intensity, sigma, snr;
- int bin, j;
-
- intensity = get_intensity(refl);
- sigma = get_esd_intensity(refl);
- snr = intensity / sigma;
-
- if ( snr < min_snr ) continue;
-
get_indices(refl, &h, &k, &l);
- res = 2.0*resolution(cell, h, k, l);
-
- bin = -1;
- for ( j=0; j<nshells; j++ ) {
- if ( (res>rmins[j]) && (res<=rmaxs[j]) ) {
- bin = j;
- break;
- }
+ if ( 2.0*resolution(cell, h, k, l) > res ) {
+ set_redundancy(refl, 0);
}
-
- /* Allow for slight rounding errors */
- if ( (bin == -1) && (res <= rmins[0]) ) bin = 0;
- if ( (bin == -1) && (res >= rmaxs[nshells-1]) ) bin = 0;
- assert(bin != -1);
-
- sh_nref[bin]++;
- sh_snr[bin] += snr;
- }
-
- for ( i=nshells-1; i>=0; i-- ) {
- if ( sh_nref[i] > 2 ) break;
}
- double cen = (rmins[i]+rmaxs[i])/2.0;
- STATUS("The highest shell is %i\n", i);
- printf("%10.3f nm^-1 or %.2f A\n", cen/1e9, 10e9/cen);
-
- FILE *fh;
- fh = fopen("cutoff.dat", "w");
- for ( i=0; i<nshells; i++ ) {
- double cen = (rmins[i]+rmaxs[i])/2.0;
- double val = sh_nref[i];
- fprintf(fh, "%10.3f %.2f %.4f\n", cen/1e9, 10e9/cen, val);
- }
- fclose(fh);
-
- free(rmins);
- free(rmaxs);
- free(sh_nref);
- free(sh_snr);
}
@@ -1819,21 +1809,28 @@ void integrate_all(struct image *image, IntegrationMethod meth,
for ( i=0; i<image->n_crystals; i++ ) {
+ double res = INFINITY;
+ Crystal *cr = image->crystals[i];
+
switch ( meth & INTEGRATION_METHOD_MASK ) {
case INTEGRATION_NONE :
break;
case INTEGRATION_RINGS :
- integrate_rings(meth, image->crystals[i], image,
+ integrate_rings(meth, cr, image,
int_diag, idh, idk, idl,
ir_inn, ir_mid, ir_out);
+ res = estimate_resolution(crystal_get_cell(cr),
+ image->features);
break;
case INTEGRATION_PROF2D :
- integrate_prof2d(meth, image->crystals[i], image,
+ integrate_prof2d(meth, cr, image,
int_diag, idh, idk, idl,
ir_inn, ir_mid, ir_out);
+ res = estimate_resolution(crystal_get_cell(cr),
+ image->features);
break;
default :
@@ -1842,10 +1839,10 @@ void integrate_all(struct image *image, IntegrationMethod meth,
}
- /* Set resolution limit */
- //resolution_cutoff(crystal_get_reflections(image->crystals[i]),
- // crystal_get_cell(image->crystals[i]));
-
+ crystal_set_resolution_limit(cr, res);
+ if ( meth & INTEGRATION_RESCUT ) {
+ apply_resolution_cutoff(cr, res);
+ }
}
}
@@ -1880,6 +1877,12 @@ IntegrationMethod integration_method(const char *str, int *err)
} else if ( strcmp(methods[i], "cen") == 0) {
meth |= INTEGRATION_CENTER;
+ } else if ( strcmp(methods[i], "rescut") == 0) {
+ meth |= INTEGRATION_RESCUT;
+
+ } else if ( strcmp(methods[i], "norescut") == 0) {
+ meth &= ~INTEGRATION_RESCUT;
+
} else if ( strcmp(methods[i], "nocen") == 0) {
meth &= ~INTEGRATION_CENTER;
diff --git a/libcrystfel/src/integration.h b/libcrystfel/src/integration.h
index e023318a..ee13b980 100644
--- a/libcrystfel/src/integration.h
+++ b/libcrystfel/src/integration.h
@@ -3,11 +3,11 @@
*
* Integration of intensities
*
- * 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 <taw@physics.org>
+ * 2010-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -56,6 +56,7 @@ typedef enum {
* @INTEGRATION_PROF2D: Two dimensional profile fitting
* @INTEGRATION_SATURATED: Integrate saturated reflections
* @INTEGRATION_CENTER: Center the peak in the box prior to integration
+ * @INTEGRATION_RESCUT: Stop integrating at the diffraction limit of the crystal
*
* An enumeration of all the available integration methods.
**/
@@ -71,6 +72,7 @@ typedef enum {
* behaviour of the integration. */
INTEGRATION_SATURATED = 256,
INTEGRATION_CENTER = 512,
+ INTEGRATION_RESCUT = 1024,
} IntegrationMethod;
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 80256829..7860209f 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -3,12 +3,12 @@
*
* Stream tools
*
- * 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 <taw@physics.org>
+ * 2010-2014 Thomas White <taw@physics.org>
* 2011 Richard Kirian
* 2011 Andrew Aquila
*
@@ -357,7 +357,7 @@ static void write_crystal(Stream *st, Crystal *cr, int include_reflections)
fprintf(st->fh, "diffraction_resolution_limit"
" = %.2f nm^-1 or %.2f A\n",
crystal_get_resolution_limit(cr)/1e9,
- 1e9 / crystal_get_resolution_limit(cr));
+ 1e10 / crystal_get_resolution_limit(cr));
fprintf(st->fh, "num_reflections = %i\n",
num_reflections(reflist));
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-17 10:07:27 +0000