aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-05-27 16:38:02 +0200
committerThomas White <taw@physics.org>2013-05-27 17:33:16 +0200
commit4c507a6c2369445092b54b2ef7859a1bc7e11803 (patch)
tree9153e8b19176ee892fdb71fd9c9b1b8349549514
parentc1c9625a8db520c685abffb886d5a43316086092 (diff)
Profile fitting stuff
-rw-r--r--libcrystfel/src/integration.c252
1 files changed, 226 insertions, 26 deletions
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index e7f61b22..e8b50930 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -150,6 +150,9 @@ struct intcontext
int n_boxes;
int max_boxes;
+ UnitCell *cell;
+ double k;
+
/* Peak region sums */
double pks_p2;
double pks_q2;
@@ -161,6 +164,7 @@ struct intcontext
int n_reference_profiles;
double **reference_profiles;
double **reference_den;
+ int *n_profiles_in_reference;
};
@@ -181,6 +185,8 @@ struct peak_box
/* Measured intensity (tentative, profile fitted or otherwise) */
double intensity;
+ double sigma;
+ double J; /* Profile scaling factor */
int rp; /* Reference profile number */
@@ -256,6 +262,8 @@ static void show_peak_box(struct intcontext *ic, struct peak_box *bx)
printf("\n");
}
printf("Reference profile number %i\n", bx->rp);
+ printf("Background parameters: a=%.2f, b=%.2f, c=%.2f\n",
+ bx->a, bx->b, bx->c);
}
@@ -264,7 +272,8 @@ static void show_reference_profile(struct intcontext *ic, int i)
{
int q;
- printf("Reference profile number %i:\n", i);
+ printf("Reference profile number %i (%i contributions):\n", i,
+ ic->n_profiles_in_reference[i]);
for ( q=ic->w-1; q>=0; q-- ) {
@@ -342,6 +351,8 @@ static void zero_profiles(struct intcontext *ic)
ic->reference_profiles[i][p+ic->w*q] = 0.0;
}
}
+
+ ic->n_profiles_in_reference[i] = 0;
}
}
@@ -433,6 +444,9 @@ static int init_intcontext(struct intcontext *ic)
if ( ic->reference_profiles == NULL ) return 1;
ic->reference_den = calloc(ic->n_reference_profiles, sizeof(double *));
if ( ic->reference_den == NULL ) return 1;
+ ic->n_profiles_in_reference = calloc(ic->n_reference_profiles,
+ sizeof(int));
+ if ( ic->n_profiles_in_reference == NULL ) return 1;
for ( i=0; i<ic->n_reference_profiles; i++ ) {
ic->reference_profiles[i] = malloc(ic->w*ic->w*sizeof(double));
if ( ic->reference_profiles[i] == NULL ) return 1;
@@ -567,6 +581,8 @@ static void add_to_reference_profile(struct intcontext *ic, struct peak_box *bx)
ic->reference_den[bx->rp][p+ic->w*q] += pow(bx->intensity, 2.0);
}
}
+
+ ic->n_profiles_in_reference[bx->rp]++;
}
@@ -605,7 +621,6 @@ static void calculate_reference_profiles(struct intcontext *ic)
}
- show_reference_profile(ic, 2);
//for ( i=0; i<ic->n_reference_profiles; i++ ) {
// show_reference_profile(ic, i);
//}
@@ -617,14 +632,27 @@ static int check_box(struct intcontext *ic, struct peak_box *bx)
int p, q;
int n_pk = 0;
int n_bg = 0;
+ double adx, ady, adz;
+ double bdx, bdy, bdz;
+ double cdx, cdy, cdz;
+ signed int hr, kr, lr;
bx->bm = malloc(ic->w*ic->w*sizeof(int));
if ( bx->bm == NULL ) return 1;
+ cell_get_cartesian(ic->cell,
+ &adx, &ady, &adz,
+ &bdx, &bdy, &bdz,
+ &cdx, &cdy, &cdz);
+ get_indices(bx->refl, &hr, &kr, &lr);
+
for ( p=0; p<ic->w; p++ ) {
for ( q=0; q<ic->w; q++ ) {
int fs, ss;
+ double hd, kd, ld;
+ signed int h, k, l;
+ struct rvec dv;
fs = bx->cfs + p;
ss = bx->css + q;
@@ -645,6 +673,19 @@ static int check_box(struct intcontext *ic, struct peak_box *bx)
bx->bm[p+ic->w*q] = BM_IG;
}
+ /* Ignore if this pixel is closer to the next reciprocal lattice
+ * point */
+ dv = get_q_for_panel(bx->p, fs, ss, NULL, ic->k);
+ hd = dv.u * adx + dv.v * ady + dv.w * adz;
+ kd = dv.u * bdx + dv.v * bdy + dv.w * bdz;
+ ld = dv.u * cdx + dv.v * cdy + dv.w * cdz;
+ h = lrint(hd);
+ k = lrint(kd);
+ l = lrint(ld);
+ if ( (h != hr) || (k != kr) || (l != lr) ) {
+ bx->bm[p+ic->w*q] = BM_IG;
+ }
+
if ( bx->bm[p+ic->w*q] == BM_PK ) n_pk++;
if ( bx->bm[p+ic->w*q] == BM_BG ) n_bg++;
@@ -658,11 +699,10 @@ static int check_box(struct intcontext *ic, struct peak_box *bx)
}
-static double fit_intensity(struct intcontext *ic, struct peak_box *bx)
+static double fit_J(struct intcontext *ic, struct peak_box *bx)
{
int p, q;
double sum = 0.0;
- double sp = 0.0;
double den = 0.0;
for ( p=0; p<ic->w; p++ ) {
@@ -677,18 +717,155 @@ static double fit_intensity(struct intcontext *ic, struct peak_box *bx)
sum += bi*P;
sum += - bx->a*p*P - bx->b*q*P - bx->c*P;
- sp += P;
den += pow(P, 2.0);
}
}
- return sum*sp / den;
+ return sum / den;
+}
+
+
+static double fit_intensity(struct intcontext *ic, struct peak_box *bx)
+{
+ int p, q;
+ double J = fit_J(ic, bx);
+ double sum = 0.0;
+
+ bx->J = J;
+
+ for ( p=0; p<ic->w; p++ ) {
+ for ( q=0; q<ic->w; q++ ) {
+
+ double P;
+
+ if ( bx->bm[p + ic->w*q] != BM_PK ) continue;
+
+ P = ic->reference_profiles[bx->rp][p+ic->w*q];
+ sum += P;
+
+ }
+ }
+
+ if ( bx->verbose ) {
+ STATUS("J = %f\n", J);
+ }
+
+ return J * sum;
}
-static void measure_all_intensities(RefList *list, struct image *image)
+static double calc_sigma(struct intcontext *ic, struct peak_box *bx)
+{
+ int p, q;
+ double sum = 0.0;
+ double mb = 0.0;
+ int nb = 0;
+ double sigb2 = 0.0;
+
+ for ( p=0; p<ic->w; p++ ) {
+ for ( q=0; q<ic->w; q++ ) {
+
+ double bi;
+
+ bi = boxi(ic, bx, p, q);
+
+ if ( bx->bm[p + ic->w*q] == BM_PK ) {
+
+ double p1, p2;
+
+ p1 = bx->J * ic->reference_profiles[bx->rp][p+ic->w*q];
+
+ p2 = bi - bx->a*p - bx->b*q - bx->c;
+ sum += pow(p1-p2, 2.0);
+
+ } else if ( bx->bm[p + ic->w*q] == BM_BG ) {
+
+ mb += bi;
+ nb++;
+
+ }
+
+ }
+ }
+
+ mb /= nb;
+
+ for ( p=0; p<ic->w; p++ ) {
+ for ( q=0; q<ic->w; q++ ) {
+
+ double bi;
+
+ if ( bx->bm[p + ic->w*q] != BM_BG ) continue;
+ bi = boxi(ic, bx, p, q);
+ sigb2 += pow(bi - mb, 2.0);
+
+ }
+ }
+
+ return sqrt(sum + sigb2);
+}
+
+
+static double average_background(struct intcontext *ic, struct peak_box *bx)
+{
+ int p, q;
+ double sum = 0.0;
+ int n = 0;
+
+ for ( p=0; p<ic->w; p++ ) {
+ for ( q=0; q<ic->w; q++ ) {
+
+ if ( bx->bm[p + ic->w*q] != BM_PK ) continue;
+
+ sum += bx->a*p + bx->b*q + bx->c;
+ n++;
+
+ }
+ }
+
+ return sum/n;
+}
+
+
+static int bg_ok(struct peak_box *bx)
+{
+ if ( (fabs(bx->a) > 10.0) || (fabs(bx->b) > 10.0) ) {
+ return 0;
+ } else {
+ return 1;
+ }
+}
+
+
+static int suitable_reference(struct intcontext *ic, struct peak_box *bx)
+{
+ int p, q;
+ double max = 0.0;
+ int height_ok;
+
+ for ( p=0; p<ic->w; p++ ) {
+ for ( q=0; q<ic->w; q++ ) {
+
+ double bi;
+
+ if ( bx->bm[p + ic->w*q] != BM_PK ) continue;
+
+ bi = boxi(ic, bx, p, q);
+ if ( bi > max ) max = bi;
+
+ }
+ }
+
+ height_ok = max > 10.0 * average_background(ic, bx);
+
+ return bg_ok(bx) && height_ok;
+}
+
+
+static void measure_all_intensities(RefList *list, struct image *image,
+ UnitCell *cell)
{
Reflection *refl;
RefListIterator *iter;
@@ -697,6 +874,8 @@ static void measure_all_intensities(RefList *list, struct image *image)
ic.halfw = 4;
ic.image = image;
+ ic.k = 1.0/image->lambda;
+ ic.cell = cell;
if ( init_intcontext(&ic) ) {
ERROR("Failed to initialise integration.\n");
return;
@@ -739,6 +918,9 @@ static void measure_all_intensities(RefList *list, struct image *image)
bx->p = p;
bx->pn = pn;
+ /* Which reference profile? */
+ bx->rp = 0;//bx->pn;
+
if ( check_box(&ic, bx) ) {
delete_box(&ic, bx);
continue;
@@ -761,9 +943,9 @@ static void measure_all_intensities(RefList *list, struct image *image)
bx->intensity = tentative_intensity(&ic, bx);
set_intensity(refl, bx->intensity);
- /* Which reference profile? */
- bx->rp = bx->pn;
- add_to_reference_profile(&ic, bx);
+ if ( suitable_reference(&ic, bx) ) {
+ add_to_reference_profile(&ic, bx);
+ }
}
calculate_reference_profiles(&ic);
@@ -774,20 +956,36 @@ static void measure_all_intensities(RefList *list, struct image *image)
bx = &ic.boxes[i];
if ( bx->verbose ) {
+ show_reference_profile(&ic, bx->rp);
STATUS("%f -> ", bx->intensity);
}
bx->intensity = fit_intensity(&ic, bx);
+ bx->sigma = calc_sigma(&ic, bx);
+
+#if 0
if ( isnan(bx->intensity) ) {
signed int h, k, l;
get_indices(bx->refl, &h, &k, &l);
- STATUS("%i %i %i\n", h, k, l);
+ STATUS("NaN intensity for %i %i %i !\n", h, k, l);
STATUS("panel %s\n", image->det->panels[bx->pn].name);
show_peak_box(&ic, bx);
+ show_reference_profile(&ic, bx->rp);
}
- set_intensity(bx->refl, bx->intensity);
- set_redundancy(bx->refl, 1);
- if ( bx->verbose ) {
- STATUS("%f\n", bx->intensity);
+ if ( bx->intensity < 0.0 ) {
+ signed int h, k, l;
+ get_indices(bx->refl, &h, &k, &l);
+ STATUS("Negative intensity (%f) for %i %i %i !\n",
+ bx->intensity, h, k, l);
+ STATUS("panel %s\n", image->det->panels[bx->pn].name);
+ show_peak_box(&ic, bx);
+ show_reference_profile(&ic, bx->rp);
+ }
+#endif
+
+ if ( bg_ok(bx) ) {
+ set_intensity(bx->refl, bx->intensity);
+ set_esd_intensity(bx->refl, bx->sigma);
+ set_redundancy(bx->refl, 1);
}
}
}
@@ -806,7 +1004,7 @@ static void estimate_mosaicity(Crystal *cr, struct image *image)
crystal_set_mosaicity(cr, mmax);
list = find_intersections(image, cr);
crystal_set_reflections(cr, list);
- measure_all_intensities(list, image);
+ measure_all_intensities(list, image, crystal_get_cell(cr));
for ( i=1; i<=msteps; i++ ) {
@@ -940,9 +1138,7 @@ static void estimate_resolution(RefList *reflections, Crystal *cr,
snr = fabs(intensity)/sigma;
/* Record intensity and set redundancy to 1 on success */
- if ( !cutoff ) {
- set_redundancy(refl, 1);
- } else {
+ if ( cutoff ) {
set_redundancy(refl, 0);
}
@@ -972,28 +1168,32 @@ static void integrate_refine(Crystal *cr, struct image *image, int use_closer,
int integrate_saturated, int **bgMasks)
{
RefList *reflections;
+ UnitCell *cell;
+
+ cell = crystal_get_cell(cr);
/* Create initial list of reflections with nominal parameters */
reflections = find_intersections(image, cr);
- measure_all_intensities(reflections, image);
+ measure_all_intensities(reflections, image, cell);
/* Find resolution limit of pattern using this list */
- estimate_resolution(reflections, cr, image);
+ //estimate_resolution(reflections, cr, image);
- reflist_free(reflections);
+ //reflist_free(reflections);
STATUS("Initial resolution estimate = %.2f nm^-1 or %.2f A\n",
crystal_get_resolution_limit(cr)/1e9,
1e9 / crystal_get_resolution_limit(cr));
/* Estimate the mosaicity of the crystal using this resolution limit */
- estimate_mosaicity(cr, image);
+ //estimate_mosaicity(cr, image);
/* Create new list of reflections with refined mosaicity */
- reflections = find_intersections(image, cr);
- measure_all_intensities(reflections, image);
+ //reflections = find_intersections(image, cr);
+ //measure_all_intensities(reflections, image);
+ crystal_set_reflections(cr, reflections);
- estimate_resolution(reflections, cr, image);
+ //estimate_resolution(reflections, cr, image);
}