diff options
| -rw-r--r-- | libcrystfel/src/integration.c | 141 |
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); } |