aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/integration.c191
1 files changed, 120 insertions, 71 deletions
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index a4cd6fc9..ff9cb7d1 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -215,6 +215,23 @@ struct intcontext
};
+struct peak_box
+{
+ int fid_fs; /* Coordinates of corner */
+ int fid_ss;
+
+ int pn; /* Panel number */
+ struct panel *p; /* The panel itself */
+
+ /* Fitted background parameters */
+ double a;
+ double b;
+ double c;
+
+ int verbose;
+};
+
+
static void addm(gsl_matrix *m, int i, int j, double val)
{
double v = gsl_matrix_get(m, i, j);
@@ -229,28 +246,68 @@ static void addv(gsl_vector *v, int i, double val)
}
-static float boxi(struct intcontext *ic, int pn,
- int fid_fs, int fid_ss, int p, int q)
+static float boxi(struct intcontext *ic, struct peak_box *bx, int p, int q)
{
int fs, ss;
- int pw;
- fs = fid_fs + p;
- ss = fid_ss + q;
+ fs = bx->fid_fs + p;
+ ss = bx->fid_ss + q;
+
+ assert(fs >= 0);
+ assert(fs < bx->p->w);
+ assert(ss >= 0);
+ assert(ss < bx->p->h);
- pw = ic->image->det->panels[pn].w;
+ assert(p >= 0);
+ assert(p < ic->w);
+ assert(q >= 0);
+ assert(q < ic->w);
- return ic->image->dp[pn][fs + pw*ss];
+ return ic->image->dp[bx->pn][fs + bx->p->w*ss];
}
-static void fit_bg(struct intcontext *ic, int pn,
- int fid_fs, int fid_ss,
- double *pa, double *pb, double *pc)
+static void show_peak_box(struct intcontext *ic, struct peak_box *bx)
+{
+ int q;
+
+ printf("Pixel values ");
+ printf("Box flags ");
+ printf("Fitted background\n");
+
+ for ( q=ic->w-1; q>=0; q-- ) {
+
+ int p;
+
+ for ( p=0; p<ic->w; p++ ) {
+ printf("%5.0f ", boxi(ic, bx, p, q));
+ }
+
+ printf(" ");
+
+ for ( p=0; p<ic->w; p++ ) {
+ printf("%i ", ic->bm[p+q*ic->w]);
+ }
+
+ printf(" ");
+
+ for ( p=0; p<ic->w; p++ ) {
+ printf("%5.0f ", bx->a*p + bx->b*q + bx->c);
+ }
+
+ printf("\n");
+ }
+ printf("-----------> p\n");
+
+}
+
+
+static void fit_bg(struct intcontext *ic, struct peak_box *bx)
{
int p, q;
gsl_vector *v;
gsl_vector *ans;
+ gsl_matrix *M;
v = gsl_vector_calloc(3);
@@ -261,7 +318,7 @@ static void fit_bg(struct intcontext *ic, int pn,
if ( ic->bm[p + ic->w*q] == BM_BG ) {
- bi = boxi(ic, pn, fid_fs, fid_ss, p, q);
+ bi = boxi(ic, bx, p, q);
addv(v, 0, bi*p);
addv(v, 1, bi*q);
@@ -272,14 +329,25 @@ static void fit_bg(struct intcontext *ic, int pn,
}
}
- ans = solve_svd(v, ic->bgm);
+ if ( bx->verbose ) {
+ show_matrix_eqn(ic->bgm, v, 3);
+ }
+
+ M = gsl_matrix_alloc(3, 3);
+ gsl_matrix_memcpy(M, ic->bgm);
+ ans = solve_svd(v, M);
gsl_vector_free(v);
+ gsl_matrix_free(M);
- *pa = gsl_vector_get(ans, 0);
- *pb = gsl_vector_get(ans, 1);
- *pc = gsl_vector_get(ans, 2);
+ bx->a = gsl_vector_get(ans, 0);
+ bx->b = gsl_vector_get(ans, 1);
+ bx->c = gsl_vector_get(ans, 2);
gsl_vector_free(ans);
+
+ if ( bx->verbose ) {
+ show_peak_box(ic, bx);
+ }
}
@@ -352,9 +420,7 @@ static int init_intcontext(struct intcontext *ic)
}
-static double tentative_intensity(struct intcontext *ic, int pn,
- int fid_fs, int fid_ss,
- double a, double b, double c)
+static double tentative_intensity(struct intcontext *ic, struct peak_box *bx)
{
int p, q;
double intensity = 0.0;
@@ -363,22 +429,20 @@ static double tentative_intensity(struct intcontext *ic, int pn,
for ( q=0; q<ic->w; q++ ) {
if ( ic->bm[p + ic->w*q] != BM_PK ) continue;
- intensity += boxi(ic, pn, fid_fs, fid_ss, p, q);
+ intensity += boxi(ic, bx, p, q);
}
}
- intensity -= a * ic->pks_p;
- intensity -= b * ic->pks_q;
- intensity -= c * ic->m;
+ intensity -= bx->a * ic->pks_p;
+ intensity -= bx->b * ic->pks_q;
+ intensity -= bx->c * ic->m;
return intensity;
}
-static void observed_position(struct intcontext *ic, int pn,
- int fid_fs, int fid_ss,
- double a, double b, double c,
+static void observed_position(struct intcontext *ic, struct peak_box *bx,
double *pos_p, double *pos_q)
{
int p, q;
@@ -391,7 +455,7 @@ static void observed_position(struct intcontext *ic, int pn,
int bi;
if ( ic->bm[p + ic->w*q] != BM_PK ) continue;
- bi = boxi(ic, pn, fid_fs, fid_ss, p, q);
+ bi = boxi(ic, bx, p, q);
num_p += bi*p;
num_q += bi*q;
@@ -400,34 +464,15 @@ static void observed_position(struct intcontext *ic, int pn,
}
}
- num_p += -a*ic->pks_p2 - b*ic->pks_pq - c*ic->pks_p;
- num_q += -a*ic->pks_q2 - b*ic->pks_pq - c*ic->pks_q;
- den += -a*ic->pks_p - b*ic->pks_q - c;
+ 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;
*pos_p = num_p / den;
*pos_q = num_q / den;
}
-static void show_peak_box(struct intcontext *ic, int pn, int fid_fs, int fid_ss)
-{
- int p, q;
-
- for ( q=ic->w-1; q>0; q-- ) {
- for ( p=0; p<ic->w; p++ ) {
- float bi;
- bi = boxi(ic, pn, fid_fs, fid_ss, p, q);
-
- printf("%5.0f ", bi);
-
- }
- printf("\n");
- }
- printf("-----------> p\n");
-
-}
-
-
static void measure_all_intensities(RefList *list, struct image *image)
{
Reflection *refl;
@@ -446,44 +491,48 @@ static void measure_all_intensities(RefList *list, struct image *image)
refl = next_refl(refl, iter) )
{
double pfs, pss;
- int fid_fs, fid_ss;
- double a, b, c;
double intensity;
- int pn;
int pw, ph;
double pos_p, pos_q;
signed int h, k, l;
+ struct peak_box bx;
+
+ bx.verbose = 0;
+
+ get_indices(refl, &h, &k, &l);
+ if ( (h==-24) && (k==6) && (l==-12) ) {
+ bx.verbose = 1;
+ }
get_detector_pos(refl, &pfs, &pss);
- fid_fs = lrint(pfs);
- fid_ss = lrint(pss);
- pn = find_panel_number(image->det, fid_fs, fid_ss);
+ bx.fid_fs = lrint(pfs);
+ bx.fid_ss = lrint(pss);
+ bx.pn = find_panel_number(image->det, bx.fid_fs, bx.fid_ss);
+ bx.p = &image->det->panels[bx.pn];
- fid_fs -= image->det->panels[pn].min_fs;
- fid_ss -= image->det->panels[pn].min_ss;
+ bx.fid_fs -= bx.p->min_fs;
+ bx.fid_ss -= bx.p->min_ss;
- fid_fs -= ic.halfw;
- fid_ss -= ic.halfw;
+ bx.fid_fs -= ic.halfw;
+ bx.fid_ss -= ic.halfw;
- pw = ic.image->det->panels[pn].w;
- ph = ic.image->det->panels[pn].h;
- if ( (fid_fs + ic.w >= pw) || (fid_ss + ic.w >= ph ) ) {
+ pw = bx.p->w;
+ ph = bx.p->h;
+ if ( (bx.fid_fs + ic.w >= pw) || (bx.fid_ss + ic.w >= ph ) ) {
+ continue;
+ }
+ if ( (bx.fid_fs < 0) || (bx.fid_ss < 0 ) ) {
continue;
}
- fit_bg(&ic, pn, fid_fs, fid_ss, &a, &b, &c);
+ fit_bg(&ic, &bx);
- get_indices(refl, &h, &k, &l);
- if ( (h==-24) && (k==6) && (l==-12) ) {
- show_peak_box(&ic, pn, fid_fs, fid_ss);
+ observed_position(&ic, &bx, &pos_p, &pos_q);
+ if ( bx.verbose ) {
+ STATUS("%f %f\n", pos_p, pos_q);
}
- observed_position(&ic, pn, fid_fs, fid_ss, a, b, c,
- &pos_p, &pos_q);
- //STATUS("%f %f\n", pos_p, pos_q);
-
- intensity = tentative_intensity(&ic, pn, fid_fs, fid_ss,
- a, b, c);
+ intensity = tentative_intensity(&ic, &bx);
set_intensity(refl, intensity);
}
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-08 17:59:07 +0000