aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2020-07-27 17:18:44 +0200
committerThomas White <taw@physics.org>2020-07-29 18:53:45 +0200
commit281b1ad167323713cc90500ce97a4224ba0272b0 (patch)
tree26b11b373b535e5eaedf231ea8631de9d5e5708f
parent55df542018664cb519e846b4db63281867be24a0 (diff)
Add low-level integration API (for better testing)
-rw-r--r--libcrystfel/src/integration.c219
-rw-r--r--libcrystfel/src/integration.h14
-rw-r--r--tests/integration_check.c27
3 files changed, 139 insertions, 121 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);
}
diff --git a/libcrystfel/src/integration.h b/libcrystfel/src/integration.h
index 629cd9ee..f402682f 100644
--- a/libcrystfel/src/integration.h
+++ b/libcrystfel/src/integration.h
@@ -109,8 +109,22 @@ typedef enum {
extern "C" {
#endif
+struct intcontext;
+
extern IntegrationMethod integration_method(const char *t, int *err);
+extern struct intcontext *intcontext_new(struct image *image,
+ UnitCell *cell,
+ IntegrationMethod meth,
+ int ir_inn, int ir_mid, int ir_out,
+ int **masks);
+
+extern int integrate_rings_once(Reflection *refl,
+ struct intcontext *ic,
+ pthread_mutex_t *term_lock);
+
+extern void intcontext_free(struct intcontext *ic);
+
extern void integrate_all(struct image *image, IntegrationMethod meth,
double ir_inn, double ir_mid, double ir_out,
IntDiag int_diag,
diff --git a/tests/integration_check.c b/tests/integration_check.c
index 8dcf5039..5fb544b7 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -31,12 +31,12 @@
#include <image.h>
#include <utils.h>
+#include <cell.h>
+#include <cell-utils.h>
+#include <integration.h>
#include "histogram.h"
-#include "../libcrystfel/src/integration.c"
-
-
int main(int argc, char *argv[])
{
struct image image;
@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
RefList *list;
Reflection *refl;
UnitCell *cell;
- struct intcontext ic;
+ struct intcontext *ic;
const int ir_inn = 2;
const int ir_mid = 4;
const int ir_out = 6;
@@ -121,25 +121,14 @@ int main(int argc, char *argv[])
deg2rad(90.0), deg2rad(90.0), deg2rad(90.0));
cell = cell_rotate(cell, random_quaternion(rng));
- 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.meth = INTEGRATION_RINGS;
- ic.int_diag = INTDIAG_NONE;
- ic.masks = NULL;
- if ( init_intcontext(&ic) ) {
+ ic = intcontext_new(&image, cell, INTEGRATION_RINGS,
+ ir_inn, ir_mid, ir_out, NULL);
+ if ( ic == NULL ) {
ERROR("Failed to initialise integration.\n");
return 1;
}
- setup_ring_masks(&ic, ir_inn, ir_mid, ir_out);
- integrate_rings_once(refl, &image, &ic, cell, 0);
+ integrate_rings_once(refl, ic, 0);
cell_free(cell);