From 281b1ad167323713cc90500ce97a4224ba0272b0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 27 Jul 2020 17:18:44 +0200 Subject: Add low-level integration API (for better testing) --- libcrystfel/src/integration.c | 219 ++++++++++++++++++++++-------------------- 1 file changed, 117 insertions(+), 102 deletions(-) (limited to 'libcrystfel/src/integration.c') diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c index 00f1d997..3c1dba0d 100644 --- a/libcrystfel/src/integration.c +++ b/libcrystfel/src/integration.c @@ -430,33 +430,106 @@ static int alloc_boxes(struct intcontext *ic, int new_max_boxes) } -static int init_intcontext(struct intcontext *ic) +static void setup_ring_masks(struct intcontext *ic, + double ir_inn, + double ir_mid, + double ir_out) { - int i; + double lim_sq, out_lim_sq, mid_lim_sq; + int p, q; + + lim_sq = pow(ir_inn, 2.0); + mid_lim_sq = pow(ir_mid, 2.0); + out_lim_sq = pow(ir_out, 2.0); + + for ( p=0; pw; p++ ) { + for ( q=0; qw; q++ ) { + + int rsq; + + rsq = (p-ic->halfw)*(p-ic->halfw) + (q-ic->halfw)*(q-ic->halfw); + + if ( rsq > out_lim_sq ) { + /* Outside outer radius */ + ic->bm[p + ic->w*q] = BM_IG; + } else { + if ( rsq >= mid_lim_sq ) { + /* Inside outer radius, outside middle radius */ + ic->bm[p + ic->w*q] = BM_BG; + } else if ( rsq <= lim_sq ) { + /* Inside inner radius */ + ic->bm[p + ic->w*q] = BM_PK; + } else { + /* Outside inner radius, inside middle radius */ + ic->bm[p + ic->w*q] = BM_IG; + } + + } + + } + } + +} + + +void intcontext_set_diag(struct intcontext *ic, + IntDiag int_diag, + signed int idh, + signed int idk, + signed int idl) +{ + ic->int_diag = int_diag; + ic->int_diag_h = idh; + ic->int_diag_k = idk; + ic->int_diag_l = idl; +} + + +struct intcontext *intcontext_new(struct image *image, + UnitCell *cell, + IntegrationMethod meth, + int ir_inn, int ir_mid, int ir_out, + int **masks) +{ + int i; + struct intcontext *ic; + + ic = malloc(sizeof(struct intcontext)); + if ( ic == NULL ) return NULL; + + ic->halfw = ir_out; + ic->image = image; + ic->k = 1.0/image->lambda; + ic->meth = meth; + ic->n_saturated = 0; + ic->n_implausible = 0; + ic->cell = cell; + ic->masks = masks; + ic->int_diag = INTDIAG_NONE; ic->w = 2*ic->halfw + 1; ic->bm = malloc(ic->w * ic->w * sizeof(enum boxmask_val)); if ( ic->bm == NULL ) { ERROR("Failed to allocate box mask.\n"); - return 1; + return NULL; } /* How many reference profiles? */ ic->n_reference_profiles = 1; ic->reference_profiles = calloc(ic->n_reference_profiles, sizeof(double *)); - if ( ic->reference_profiles == NULL ) return 1; + if ( ic->reference_profiles == NULL ) return NULL; ic->reference_den = calloc(ic->n_reference_profiles, sizeof(double *)); - if ( ic->reference_den == NULL ) return 1; + if ( ic->reference_den == NULL ) return NULL; ic->n_profiles_in_reference = calloc(ic->n_reference_profiles, sizeof(int)); - if ( ic->n_profiles_in_reference == NULL ) return 1; + if ( ic->n_profiles_in_reference == NULL ) return NULL; for ( i=0; in_reference_profiles; i++ ) { ic->reference_profiles[i] = malloc(ic->w*ic->w*sizeof(double)); - if ( ic->reference_profiles[i] == NULL ) return 1; + if ( ic->reference_profiles[i] == NULL ) return NULL; ic->reference_den[i] = malloc(ic->w*ic->w*sizeof(double)); - if ( ic->reference_den[i] == NULL ) return 1; + if ( ic->reference_den[i] == NULL ) return NULL; } zero_profiles(ic); @@ -464,14 +537,16 @@ static int init_intcontext(struct intcontext *ic) ic->n_boxes = 0; ic->max_boxes = 0; if ( alloc_boxes(ic, 32) ) { - return 1; + return NULL; } - return 0; + setup_ring_masks(ic, ir_inn, ir_mid, ir_out); + + return ic; } -static void free_intcontext(struct intcontext *ic) +void intcontext_free(struct intcontext *ic) { int i; @@ -492,47 +567,6 @@ static void free_intcontext(struct intcontext *ic) } -static void setup_ring_masks(struct intcontext *ic, - double ir_inn, double ir_mid, double ir_out) -{ - double lim_sq, out_lim_sq, mid_lim_sq; - int p, q; - - lim_sq = pow(ir_inn, 2.0); - mid_lim_sq = pow(ir_mid, 2.0); - out_lim_sq = pow(ir_out, 2.0); - - for ( p=0; pw; p++ ) { - for ( q=0; qw; q++ ) { - - int rsq; - - rsq = (p-ic->halfw)*(p-ic->halfw) + (q-ic->halfw)*(q-ic->halfw); - - if ( rsq > out_lim_sq ) { - /* Outside outer radius */ - ic->bm[p + ic->w*q] = BM_IG; - } else { - - if ( rsq >= mid_lim_sq ) { - /* Inside outer radius, outside middle radius */ - ic->bm[p + ic->w*q] = BM_BG; - } else if ( rsq <= lim_sq ) { - /* Inside inner radius */ - ic->bm[p + ic->w*q] = BM_PK; - } else { - /* Outside inner radius, inside middle radius */ - ic->bm[p + ic->w*q] = BM_IG; - } - - } - - } - } - -} - - static struct peak_box *add_box(struct intcontext *ic) { int idx; @@ -1346,55 +1380,46 @@ static void integrate_prof2d(IntegrationMethod meth, { RefList *list; UnitCell *cell; - struct intcontext ic; + struct intcontext *ic; int i; list = crystal_get_reflections(cr); cell = crystal_get_cell(cr); - ic.halfw = ir_out; - ic.image = image; - ic.k = 1.0/image->lambda; - ic.meth = meth; - ic.n_saturated = 0; - ic.n_implausible = 0; - ic.cell = cell; - ic.int_diag = int_diag; - ic.int_diag_h = idh; - ic.int_diag_k = idk; - ic.int_diag_l = idl; - ic.masks = masks; - if ( init_intcontext(&ic) ) { + ic = intcontext_new(image, cell, meth, + ir_inn, ir_mid, ir_out, + masks); + if ( ic == NULL ) { ERROR("Failed to initialise integration.\n"); return; } - setup_ring_masks(&ic, ir_inn, ir_mid, ir_out); - setup_profile_boxes(&ic, list); - calculate_reference_profiles(&ic); + intcontext_set_diag(ic, int_diag, idh, idk, idl); + setup_profile_boxes(ic, list); + calculate_reference_profiles(ic); - for ( i=0; in_reference_profiles; i++ ) { + if ( ic->n_profiles_in_reference[i] == 0 ) { ERROR("Reference profile %i has no contributions.\n", i); - free_intcontext(&ic); + intcontext_free(ic); return; } } - for ( i=0; in_boxes; i++ ) { struct peak_box *bx; - bx = &ic.boxes[i]; - integrate_prof2d_once(&ic, bx, term_lock); + bx = &ic->boxes[i]; + integrate_prof2d_once(ic, bx, term_lock); } - free_intcontext(&ic); + intcontext_free(ic); } -static int integrate_rings_once(Reflection *refl, struct image *image, - struct intcontext *ic, UnitCell *cell, - pthread_mutex_t *term_lock) +int integrate_rings_once(Reflection *refl, + struct intcontext *ic, + pthread_mutex_t *term_lock) { double pfs, pss; struct peak_box *bx; @@ -1429,7 +1454,7 @@ static int integrate_rings_once(Reflection *refl, struct image *image, bx->refl = refl; bx->cfs = cfs; bx->css = css; - bx->p = &image->detgeom->panels[pn]; + bx->p = &ic->image->detgeom->panels[pn]; bx->pn = pn; if ( ic->meth & INTEGRATION_CENTER ) { @@ -1606,49 +1631,39 @@ static void integrate_rings(IntegrationMethod meth, Reflection *refl; RefListIterator *iter; UnitCell *cell; - struct intcontext ic; + struct intcontext *ic; int n_rej = 0; list = crystal_get_reflections(cr); cell = crystal_get_cell(cr); - ic.halfw = ir_out; - ic.image = image; - ic.k = 1.0/image->lambda; - ic.n_saturated = 0; - ic.n_implausible = 0; - ic.cell = cell; - ic.ir_inn = ir_inn; - ic.ir_mid = ir_mid; - ic.ir_out = ir_out; - ic.int_diag = int_diag; - ic.int_diag_h = idh; - ic.int_diag_k = idk; - ic.int_diag_l = idl; - ic.meth = meth; - ic.masks = masks; - if ( init_intcontext(&ic) ) { + ic = intcontext_new(image, cell, meth, + ir_inn, ir_mid, ir_out, + masks); + if ( ic == NULL ) { ERROR("Failed to initialise integration.\n"); return; } - setup_ring_masks(&ic, ir_inn, ir_mid, ir_out); + + intcontext_set_diag(ic, int_diag, idh, idk, idl); for ( refl = first_refl(list, &iter); refl != NULL; refl = next_refl(refl, iter) ) { - n_rej += integrate_rings_once(refl, image, &ic, cell, term_lock); + n_rej += integrate_rings_once(refl, ic, + term_lock); } - free_intcontext(&ic); + intcontext_free(ic); if ( n_rej > 0 ) { ERROR("WARNING: %i reflections could not be integrated\n", n_rej); } - crystal_set_num_saturated_reflections(cr, ic.n_saturated); - crystal_set_num_implausible_reflections(cr, ic.n_implausible); + crystal_set_num_saturated_reflections(cr, ic->n_saturated); + crystal_set_num_implausible_reflections(cr, ic->n_implausible); } -- cgit v1.2.3