aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/integration.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-06-21 17:35:20 +0200
committerThomas White <taw@physics.org>2013-06-21 17:35:20 +0200
commit5683a03c1aecb02688fe95145bfa7670577ce00d (patch)
tree23d021c1fb971dc9c1591235a27a5e7ccaa1e1a2 /libcrystfel/src/integration.c
parentc0b30bb042b8c5304efe31f6af0baf827409c4f1 (diff)
Peak integrals are per box
Diffstat (limited to 'libcrystfel/src/integration.c')
-rw-r--r--libcrystfel/src/integration.c148
1 files changed, 75 insertions, 73 deletions
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index 3cf2ff97..f00e2227 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -156,7 +156,6 @@ struct intcontext
int w;
enum boxmask_val *bm; /* Box mask */
struct image *image;
- gsl_matrix *bgm; /* Background estimation matrix */
struct peak_box *boxes;
int n_boxes;
@@ -165,14 +164,6 @@ struct intcontext
UnitCell *cell;
double k;
- /* Peak region sums */
- double pks_p2;
- double pks_q2;
- double pks_pq;
- double pks_p;
- double pks_q;
- int m;
-
int n_reference_profiles;
double **reference_profiles;
double **reference_den;
@@ -199,6 +190,16 @@ struct peak_box
double b;
double c;
+ /* Peak region sums */
+ double pks_p2;
+ double pks_q2;
+ double pks_pq;
+ double pks_p;
+ double pks_q;
+ int m;
+
+ gsl_matrix *bgm; /* Background estimation matrix */
+
/* Measured intensity (tentative, profile fitted or otherwise) */
double intensity;
double sigma;
@@ -410,11 +411,11 @@ static void fit_bg(struct intcontext *ic, struct peak_box *bx)
}
if ( bx->verbose ) {
- show_matrix_eqn(ic->bgm, v, 3);
+ show_matrix_eqn(bx->bgm, v, 3);
}
/* SVD is massive overkill here */
- ans = solve_svd(v, ic->bgm);
+ ans = solve_svd(v, bx->bgm);
gsl_vector_free(v);
bx->a = gsl_vector_get(ans, 0);
@@ -464,7 +465,6 @@ static int alloc_boxes(struct intcontext *ic, int new_max_boxes)
static int init_intcontext(struct intcontext *ic)
{
- int p, q;
int i;
ic->w = 2*ic->halfw + 1;
@@ -475,60 +475,6 @@ static int init_intcontext(struct intcontext *ic)
return 1;
}
- ic->bgm = gsl_matrix_calloc(3, 3);
- if ( ic->bgm == NULL ) {
- ERROR("Failed to initialise matrix.\n");
- return 1;
- }
-
- ic->pks_p2 = 0.0;
- ic->pks_q2 = 0.0;
- ic->pks_pq = 0.0;
- ic->pks_p = 0.0;
- ic->pks_q = 0.0;
- ic->m = 0;
-
- for ( p=0; p<ic->w; p++ ) {
- for ( q=0; q<ic->w; q++ ) {
-
- if ( (p==0) || (q==0) || (p==ic->w-1) || (q==ic->w-1) ) {
- ic->bm[p + ic->w*q] = BM_BG;
- } else {
- ic->bm[p + ic->w*q] = BM_PK;
- }
-
- switch ( ic->bm[p + ic->w*q] ) {
-
- case BM_IG :
- case BM_BH :
- break;
-
- case BM_BG :
- addm(ic->bgm, 0, 0, p*p);
- addm(ic->bgm, 0, 1, p*q);
- addm(ic->bgm, 0, 2, p);
- addm(ic->bgm, 1, 0, p*q);
- addm(ic->bgm, 1, 1, q*q);
- addm(ic->bgm, 1, 2, q);
- addm(ic->bgm, 2, 0, p);
- addm(ic->bgm, 2, 1, q);
- addm(ic->bgm, 2, 2, 1);
- break;
-
- case BM_PK :
- ic->pks_p2 += p*p;
- ic->pks_q2 += q*q;
- ic->pks_pq += p*q;
- ic->pks_p += p;
- ic->pks_q += q;
- ic->m++;
- break;
-
- }
-
- }
- }
-
/* How many reference profiles? */
ic->n_reference_profiles = ic->image->det->n_panels;
ic->reference_profiles = calloc(ic->n_reference_profiles,
@@ -564,6 +510,7 @@ static void free_intcontext(struct intcontext *ic)
for ( i=0; i<ic->n_boxes; i++ ) {
free(ic->boxes[i].bm);
+ gsl_matrix_free(ic->boxes[i].bgm);
}
free(ic->boxes);
@@ -575,7 +522,6 @@ static void free_intcontext(struct intcontext *ic)
free(ic->reference_den);
free(ic->n_profiles_in_reference);
free(ic->bm);
- gsl_matrix_free(ic->bgm);
}
@@ -647,6 +593,12 @@ static struct peak_box *add_box(struct intcontext *ic)
ic->boxes[idx].refl = NULL;
ic->boxes[idx].verbose = 0;
+ ic->boxes[idx].bgm = gsl_matrix_calloc(3, 3);
+ if ( ic->boxes[idx].bgm == NULL ) {
+ ERROR("Failed to initialise matrix.\n");
+ return NULL;
+ }
+
return &ic->boxes[idx];
}
@@ -690,9 +642,9 @@ static double tentative_intensity(struct intcontext *ic, struct peak_box *bx)
}
}
- intensity -= bx->a * ic->pks_p;
- intensity -= bx->b * ic->pks_q;
- intensity -= bx->c * ic->m;
+ intensity -= bx->a * bx->pks_p;
+ intensity -= bx->b * bx->pks_q;
+ intensity -= bx->c * bx->m;
return intensity;
}
@@ -720,9 +672,9 @@ static void observed_position(struct intcontext *ic, struct peak_box *bx,
}
}
- num_p += -bx->a*ic->pks_p2 - bx->b*ic->pks_pq - bx->c*ic->pks_p;
- num_q += -bx->a*ic->pks_q2 - bx->b*ic->pks_pq - bx->c*ic->pks_q;
- den += -bx->a*ic->pks_p - bx->b*ic->pks_q - bx->c;
+ num_p += -bx->a*bx->pks_p2 - bx->b*bx->pks_pq - bx->c*bx->pks_p;
+ num_q += -bx->a*bx->pks_q2 - bx->b*bx->pks_pq - bx->c*bx->pks_q;
+ den += -bx->a*bx->pks_p - bx->b*bx->pks_q - bx->c;
*pos_p = num_p / den;
*pos_q = num_q / den;
@@ -796,6 +748,54 @@ static void calculate_reference_profiles(struct intcontext *ic)
}
+static void setup_peak_integrals(struct intcontext *ic, struct peak_box *bx)
+{
+ int p, q;
+
+ bx->pks_p2 = 0.0;
+ bx->pks_q2 = 0.0;
+ bx->pks_pq = 0.0;
+ bx->pks_p = 0.0;
+ bx->pks_q = 0.0;
+ bx->m = 0;
+
+ for ( p=0; p<ic->w; p++ ) {
+ for ( q=0; q<ic->w; q++ ) {
+
+ switch ( bx->bm[p + ic->w*q] ) {
+
+ case BM_IG :
+ case BM_BH :
+ break;
+
+ case BM_BG :
+ addm(bx->bgm, 0, 0, p*p);
+ addm(bx->bgm, 0, 1, p*q);
+ addm(bx->bgm, 0, 2, p);
+ addm(bx->bgm, 1, 0, p*q);
+ addm(bx->bgm, 1, 1, q*q);
+ addm(bx->bgm, 1, 2, q);
+ addm(bx->bgm, 2, 0, p);
+ addm(bx->bgm, 2, 1, q);
+ addm(bx->bgm, 2, 2, 1);
+ break;
+
+ case BM_PK :
+ bx->pks_p2 += p*p;
+ bx->pks_q2 += q*q;
+ bx->pks_pq += p*q;
+ bx->pks_p += p;
+ bx->pks_q += q;
+ bx->m++;
+ break;
+
+ }
+
+ }
+ }
+}
+
+
static int check_box(struct intcontext *ic, struct peak_box *bx, int *sat)
{
int p, q;
@@ -869,6 +869,8 @@ static int check_box(struct intcontext *ic, struct peak_box *bx, int *sat)
if ( n_pk < 4 ) return 1;
if ( n_bg < 4 ) return 1;
+ setup_peak_integrals(ic, bx);
+
return 0;
}