aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-05-24 11:05:32 +0200
committerThomas White <taw@physics.org>2013-05-27 17:33:15 +0200
commitb990cd2c63a68433ed780f1f3400601bec759c5b (patch)
tree70b095885a0dfff003fb6576c7b99fc8d86ae837
parentebace165026c40452c5fdd2a8927c12b3ded5e4f (diff)
Calculate the reference profiles
-rw-r--r--libcrystfel/src/integration.c141
1 files changed, 138 insertions, 3 deletions
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index ff9cb7d1..f4cab6eb 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -212,6 +212,10 @@ struct intcontext
double pks_p;
double pks_q;
int m;
+
+ int n_reference_profiles;
+ double **reference_profiles;
+ double **reference_den;
};
@@ -228,6 +232,11 @@ struct peak_box
double b;
double c;
+ /* Measured intensity (tentative, profile fitted or otherwise) */
+ double intensity;
+
+ int rp; /* Reference profile number */
+
int verbose;
};
@@ -298,6 +307,29 @@ static void show_peak_box(struct intcontext *ic, struct peak_box *bx)
printf("\n");
}
printf("-----------> p\n");
+ printf("Reference profile number %i\n", bx->rp);
+
+}
+
+
+static void show_reference_profile(struct intcontext *ic, int i)
+{
+ int q;
+
+ printf("Reference profile number %i:\n", i);
+ printf("Pixel values\n");
+
+ for ( q=ic->w-1; q>=0; q-- ) {
+
+ int p;
+
+ for ( p=0; p<ic->w; p++ ) {
+ printf("%5.0f ", ic->reference_profiles[i][p+ic->w*q]);
+ }
+
+ printf("\n");
+ }
+ printf("-----------> p\n");
}
@@ -351,9 +383,27 @@ static void fit_bg(struct intcontext *ic, struct peak_box *bx)
}
+static void zero_profiles(struct intcontext *ic)
+{
+ int i;
+
+ for ( i=0; i<ic->n_reference_profiles; i++ ) {
+
+ int p, q;
+
+ for ( p=0; p<ic->w; p++ ) {
+ for ( q=0; q<ic->w; q++ ) {
+ ic->reference_profiles[i][p+ic->w*q] = 0.0;
+ }
+ }
+ }
+}
+
+
static int init_intcontext(struct intcontext *ic)
{
int p, q;
+ int i;
ic->w = 2*ic->halfw + 1;
@@ -416,6 +466,21 @@ static int init_intcontext(struct intcontext *ic)
}
}
+ /* How many reference profiles? */
+ ic->n_reference_profiles = ic->image->det->n_panels;
+ ic->reference_profiles = calloc(ic->n_reference_profiles,
+ sizeof(double *));
+ if ( ic->reference_profiles == NULL ) return 1;
+ ic->reference_den = calloc(ic->n_reference_profiles, sizeof(double *));
+ if ( ic->reference_den == 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;
+ ic->reference_den[i] = malloc(ic->w*ic->w*sizeof(double));
+ if ( ic->reference_den[i] == NULL ) return 1;
+ }
+ zero_profiles(ic);
+
return 0;
}
@@ -473,6 +538,69 @@ static void observed_position(struct intcontext *ic, struct peak_box *bx,
}
+static void add_to_reference_profile(struct intcontext *ic, struct peak_box *bx)
+{
+ int p, q;
+
+ for ( p=0; p<ic->w; p++ ) {
+ for ( q=0; q<ic->w; q++ ) {
+
+ double val;
+ float bi;
+
+ if ( ic->bm[p + ic->w*q] == BM_IG ) continue;
+ bi = boxi(ic, bx, p, q);
+
+ val = bi*bx->intensity;
+ val -= p*bx->intensity*bx->a + q*bx->intensity*bx->b
+ + bx->intensity*bx->c;
+
+ ic->reference_profiles[bx->rp][p+ic->w*q] += val;
+ ic->reference_den[bx->rp][p+ic->w*q] += pow(bx->intensity, 2.0);
+ }
+ }
+}
+
+
+static void calculate_reference_profiles(struct intcontext *ic)
+{
+ int i;
+
+ for ( i=0; i<ic->n_reference_profiles; i++ ) {
+
+ int p, q;
+ double max = 0.0;
+
+ for ( p=0; p<ic->w; p++ ) {
+ for ( q=0; q<ic->w; q++ ) {
+
+ double den;
+
+ den = ic->reference_den[i][p+ic->w*q];
+ ic->reference_profiles[i][p+ic->w*q] /= den;
+ if ( ic->reference_profiles[i][p+ic->w*q] > max ) {
+ max = ic->reference_profiles[i][p+ic->w*q];
+ }
+
+ }
+ }
+
+ max /= 10000.0;
+
+ for ( p=0; p<ic->w; p++ ) {
+ for ( q=0; q<ic->w; q++ ) {
+
+ ic->reference_profiles[i][p+ic->w*q] /= max;
+
+ }
+ }
+
+ }
+
+ show_reference_profile(ic, 2);
+}
+
+
static void measure_all_intensities(RefList *list, struct image *image)
{
Reflection *refl;
@@ -491,7 +619,6 @@ static void measure_all_intensities(RefList *list, struct image *image)
refl = next_refl(refl, iter) )
{
double pfs, pss;
- double intensity;
int pw, ph;
double pos_p, pos_q;
signed int h, k, l;
@@ -528,13 +655,21 @@ static void measure_all_intensities(RefList *list, struct image *image)
fit_bg(&ic, &bx);
observed_position(&ic, &bx, &pos_p, &pos_q);
+ pos_p -= ic.halfw;
+ pos_q -= ic.halfw;
if ( bx.verbose ) {
STATUS("%f %f\n", pos_p, pos_q);
}
- intensity = tentative_intensity(&ic, &bx);
- set_intensity(refl, intensity);
+ bx.intensity = tentative_intensity(&ic, &bx);
+ set_intensity(refl, bx.intensity);
+
+ /* Which reference profile? */
+ bx.rp = bx.pn;
+ add_to_reference_profile(&ic, &bx);
}
+
+ calculate_reference_profiles(&ic);
}