Add low-level integration API (for better testing)

1 files changed, 117 insertions, 102 deletions

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; p<ic->w; p++ ) {
+	for ( q=0; q<ic->w; 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; 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;
+		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; p<ic->w; p++ ) {
-	for ( q=0; q<ic->w; 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; i<ic.n_reference_profiles; i++ ) {
-		if ( ic.n_profiles_in_reference[i] == 0 ) {
+	for ( i=0; i<ic->n_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; i<ic.n_boxes; i++ ) {
+	for ( i=0; i<ic->n_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);
 }