aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/geometry.c254
-rw-r--r--libcrystfel/src/geometry.h10
-rw-r--r--libcrystfel/src/integration.c4
-rw-r--r--libcrystfel/src/reflist-utils.c2
-rw-r--r--libcrystfel/src/reflist.c45
-rw-r--r--libcrystfel/src/reflist.h10
-rw-r--r--libcrystfel/src/stream.c10
7 files changed, 129 insertions, 206 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index b55b696d..35492fde 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -47,58 +47,66 @@
#include "geometry.h"
-static signed int locate_peak(double x, double y, double z, double k,
- struct detector *det, double *xdap, double *ydap)
+static int locate_peak_on_panel(double x, double y, double z, double k,
+ struct panel *p,
+ double *pfs, double *pss)
{
- int i;
- signed int found = -1;
const double den = k + z;
+ double fs, ss, plx, ply, xd, yd;
- *xdap = -1; *ydap = -1;
+ /* Coordinates of peak relative to central beam, in m */
+ xd = p->clen * x / den;
+ yd = p->clen * y / den;
- for ( i=0; i<det->n_panels; i++ ) {
+ /* Convert to pixels */
+ xd *= p->res;
+ yd *= p->res;
- double xd, yd;
- double fs, ss, plx, ply;
- struct panel *p;
+ /* Convert to relative to the panel corner */
+ plx = xd - p->cnx;
+ ply = yd - p->cny;
+
+ fs = p->xfs*plx + p->yfs*ply;
+ ss = p->xss*plx + p->yss*ply;
+
+ fs += p->min_fs;
+ ss += p->min_ss;
+
+ *pfs = fs; *pss = ss;
+
+ /* Now, is this on this panel? */
+ if ( fs < p->min_fs ) return 0;
+ if ( fs > p->max_fs ) return 0;
+ if ( ss < p->min_ss ) return 0;
+ if ( ss > p->max_ss ) return 0;
- p = &det->panels[i];
- /* Coordinates of peak relative to central beam, in m */
- xd = p->clen * x / den;
- yd = p->clen * y / den;
+ return 1;
+}
- /* Convert to pixels */
- xd *= p->res;
- yd *= p->res;
+static signed int locate_peak(double x, double y, double z, double k,
+ struct detector *det, double *pfs, double *pss)
+{
+ int i;
- /* Convert to relative to the panel corner */
- plx = xd - p->cnx;
- ply = yd - p->cny;
+ *pfs = -1; *pss = -1;
- fs = p->xfs*plx + p->yfs*ply;
- ss = p->xss*plx + p->yss*ply;
+ for ( i=0; i<det->n_panels; i++ ) {
- fs += p->min_fs;
- ss += p->min_ss;
+ struct panel *p;
+
+ p = &det->panels[i];
- /* Now, is this on this panel? */
- if ( fs < p->min_fs ) continue;
- if ( fs > p->max_fs ) continue;
- if ( ss < p->min_ss ) continue;
- if ( ss > p->max_ss ) continue;
+ if ( locate_peak_on_panel(x, y, z, k, p, pfs, pss) ) {
- /* If peak appears on multiple panels, reject it */
- if ( found != -1 ) return -1;
+ /* Woohoo! */
+ return i;
- /* Woohoo! */
- found = i;
- *xdap = fs;
- *ydap = ss;
+ }
}
- return found;
+ return -1;
}
@@ -182,7 +190,8 @@ static double partiality(PartialityModel pmodel,
static Reflection *check_reflection(struct image *image, Crystal *cryst,
PartialityModel pmodel,
signed int h, signed int k, signed int l,
- double xl, double yl, double zl)
+ double xl, double yl, double zl,
+ Reflection *updateme)
{
const int output = 0;
double tl;
@@ -195,7 +204,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
double del;
/* Don't predict 000 */
- if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL;
+ if ( (updateme == NULL) && (abs(h)+abs(k)+abs(l) == 0) ) return NULL;
pr = crystal_get_profile_radius(cryst);
del = image->div + crystal_get_mosaicity(cryst);
@@ -207,7 +216,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0);
/* If the point is looking "backscattery", reject it straight away */
- if ( zl < -khigh/2.0 ) return NULL;
+ if ( (updateme == NULL) && (zl < -khigh/2.0) ) return NULL;
tl = sqrt(xl*xl + yl*yl);
@@ -220,28 +229,47 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */
/* Condition for reflection to be excited at all */
- if ( (signbit(rlow) == signbit(rhigh))
+ if ( (updateme == NULL)
+ && (signbit(rlow) == signbit(rhigh))
&& (fabs(rlow) > pr)
&& (fabs(rhigh) > pr) ) return NULL;
/* Calculate partiality */
part = partiality(pmodel, rlow, rhigh, pr);
- /* Add peak to list */
- refl = reflection_new(h, k, l);
+ if ( updateme == NULL ) {
+ refl = reflection_new(h, k, l);
+ } else {
+ refl = updateme;
+ }
+
+ /* If we are updating a previous reflection, assume it stays
+ * on the same panel and calculate the new position even if it's
+ * fallen off the edge of the panel. */
+ if ( updateme != NULL ) {
+
+ double fs, ss;
+ locate_peak_on_panel(xl, yl, zl, 1.0/image->lambda,
+ get_panel(updateme), &fs, &ss);
+ set_detector_pos(refl, fs, ss);
+
+ }
- /* If we have detector information, check the spot is measured.
- * Otherwise, we make do with calculating the partialiaty etc. */
- if ( image->det != NULL ) {
- double xda, yda; /* Position on detector */
- signed int p; /* Panel number */
- p = locate_peak(xl, yl, zl, 1.0/image->lambda, image->det,
- &xda, &yda);
+ /* Otherwise, calculate position if we have a detector structure, and
+ * if we don't then just make do with partiality calculation */
+ if ( (image->det != NULL) && (updateme == NULL) ) {
+
+ double fs, ss; /* Position on detector */
+ signed int p; /* Panel number */
+ p = locate_peak(xl, yl, zl, 1.0/image->lambda,
+ image->det, &fs, &ss);
if ( p == -1 ) {
reflection_free(refl);
return NULL;
}
- set_detector_pos(refl, 0.0, xda, yda);
+ set_detector_pos(refl, fs, ss);
+ set_panel(refl, &image->det->panels[p]);
+
}
if ( unlikely(rlow < rhigh) ) {
@@ -249,8 +277,11 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n",
h, k, l, rlow, rhigh);
ERROR("div + m = %e, R = %e, bw = %e\n", del, pr, image->bw);
- reflection_free(refl);
- return NULL;
+ /* If we are updating, this is (kind of) OK */
+ if ( updateme == NULL ) {
+ reflection_free(refl);
+ return NULL;
+ }
}
set_partial(refl, rlow, rhigh, part);
@@ -339,7 +370,7 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
zl = h*asz + k*bsz + l*csz;
refl = check_reflection(image, cryst, pmodel,
- h, k, l, xl, yl, zl);
+ h, k, l, xl, yl, zl, NULL);
if ( refl != NULL ) {
add_refl_to_list(refl, reflections);
@@ -353,67 +384,6 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
}
-/* Deprecated: select reflections using Kirian-style pixel proximity */
-RefList *select_intersections(struct image *image, Crystal *cryst)
-{
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- const double min_dist = 0.25;
- RefList *list;
- int i;
-
- /* Round towards nearest */
- fesetround(1);
-
- /* Cell basis vectors for this image */
- cell_get_cartesian(crystal_get_cell(cryst), &ax, &ay, &az,
- &bx, &by, &bz, &cx, &cy, &cz);
-
- list = reflist_new();
- if ( list == NULL ) return NULL;
-
- /* Loop over peaks, checking proximity to nearest reflection */
- for ( i=0; i<image_feature_count(image->features); i++ ) {
-
- struct imagefeature *f;
- struct rvec q;
- double h, k, l, hd, kd, ld;
- double dsq;
-
- f = image_get_feature(image->features, i);
- if ( f == NULL ) continue;
-
- /* Reciprocal space position of found peak */
- q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
-
- /* Decimal and fractional Miller indices of nearest
- * reciprocal lattice point */
- hd = q.u * ax + q.v * ay + q.w * az;
- kd = q.u * bx + q.v * by + q.w * bz;
- ld = q.u * cx + q.v * cy + q.w * cz;
- h = lrint(hd);
- k = lrint(kd);
- l = lrint(ld);
-
- /* Check distance */
- dsq = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
-
- if ( sqrt(dsq) < min_dist ) {
-
- Reflection *refl;
-
- refl = add_refl(list, h, k, l);
- set_detector_pos(refl, sqrt(dsq), f->fs, f->ss);
-
- }
-
- }
-
- return list;
-}
-
-
static void set_unity_partialities(Crystal *cryst)
{
Reflection *refl;
@@ -430,8 +400,7 @@ static void set_unity_partialities(Crystal *cryst)
/* Calculate partialities and apply them to the image's reflections */
-void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
- int *n_gained, int *n_lost, double *mean_p_change)
+void update_partialities(Crystal *cryst, PartialityModel pmodel)
{
Reflection *refl;
RefListIterator *iter;
@@ -439,8 +408,6 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
double bsx, bsy, bsz;
double csx, csy, csz;
struct image *image = crystal_get_image(cryst);
- double total_p_change = 0.0;
- int n = 0;
if ( pmodel == PMODEL_UNITY ) {
set_unity_partialities(cryst);
@@ -454,69 +421,20 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
refl != NULL;
refl = next_refl(refl, iter) )
{
- Reflection *vals;
- double r1, r2, L, p, x, y;
double xl, yl, zl;
signed int h, k, l;
- double old_p;
get_symmetric_indices(refl, &h, &k, &l);
- old_p = get_partiality(refl);
/* Get the coordinates of the reciprocal lattice point */
xl = h*asx + k*bsx + l*csx;
yl = h*asy + k*bsy + l*csy;
zl = h*asz + k*bsz + l*csz;
- vals = check_reflection(image, cryst, pmodel,
- h, k, l, xl, yl, zl);
-
- if ( vals == NULL ) {
-
- if ( get_redundancy(refl) != 0 ) {
- (*n_lost)++;
- set_partiality(refl, 0.0);
- set_redundancy(refl, 0);
- }
-
- } else {
-
- if ( get_redundancy(refl) == 0 ) {
- (*n_gained)++;
- set_redundancy(refl, 1);
- }
-
- /* Transfer partiality stuff */
- get_partial(vals, &r1, &r2, &p);
- set_partial(refl, r1, r2, p);
- L = get_lorentz(vals);
- set_lorentz(refl, L);
-
- /* Transfer detector location */
- get_detector_pos(vals, &x, &y);
- set_detector_pos(refl, 0.0, x, y);
-
- reflection_free(vals);
-
- total_p_change += fabs(p - old_p);
- n++;
-
- }
+ check_reflection(image, cryst, pmodel,
+ h, k, l, xl, yl, zl, refl);
}
-
- *mean_p_change = total_p_change / n;
-}
-
-
-/* Wrapper to maintain API compatibility */
-void update_partialities(Crystal *cryst, PartialityModel pmodel)
-{
- int n_gained = 0;
- int n_lost = 0;
- double mean_p_change = 0.0;
- update_partialities_2(cryst, pmodel, &n_gained, &n_lost,
- &mean_p_change);
}
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index f39b6f8b..162c0355 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -3,12 +3,12 @@
*
* Geometry of diffraction
*
- * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
* Authors:
- * 2010-2014 Thomas White <taw@physics.org>
+ * 2010-2015 Thomas White <taw@physics.org>
* 2012 Richard Kirian
*
* This file is part of CrystFEL.
@@ -67,13 +67,7 @@ extern RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
PartialityModel pmodel,
double max_res);
-/* Deprecated: select reflections using Kirian-style pixel proximity */
-extern RefList *select_intersections(struct image *image, Crystal *cryst);
-
extern void update_partialities(Crystal *cryst, PartialityModel pmodel);
-extern void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
- int *n_gained, int *n_lost,
- double *mean_p_change);
extern void polarisation_correction(RefList *list, UnitCell *cell,
struct image *image);
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index a05f64f1..0e35037d 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -1193,7 +1193,7 @@ static void integrate_prof2d_once(struct intcontext *ic, struct peak_box *bx,
get_detector_pos(bx->refl, &pfs, &pss);
pfs += bx->offs_fs;
pss += bx->offs_ss;
- set_detector_pos(bx->refl, 0.0, pfs, pss);
+ set_detector_pos(bx->refl, pfs, pss);
if ( bx->intensity < -5.0*bx->sigma ) {
ic->n_implausible++;
@@ -1447,7 +1447,7 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
/* Update position */
pfs += bx->offs_fs;
pss += bx->offs_ss;
- set_detector_pos(refl, 0.0, pfs, pss);
+ set_detector_pos(refl, pfs, pss);
if ( get_int_diag(ic, refl) ) show_peak_box(ic, bx, results_pipe);
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
index 722f800c..82943158 100644
--- a/libcrystfel/src/reflist-utils.c
+++ b/libcrystfel/src/reflist-utils.c
@@ -390,7 +390,7 @@ RefList *read_reflections_from_file(FILE *fh)
refl = add_refl(out, h, k, l);
set_intensity(refl, intensity);
- set_detector_pos(refl, 0.0, fs, ss);
+ set_detector_pos(refl, fs, ss);
set_esd_intensity(refl, sigma);
set_redundancy(refl, cts);
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index 33c1f948..795894a1 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -75,10 +75,7 @@ struct _refldata {
/* Location in image */
double fs;
double ss;
-
- /* The distance from the exact Bragg position to the coordinates
- * given above. */
- double excitation_error;
+ struct panel *panel;
/* Non-zero if this reflection can be used for scaling */
int scalable;
@@ -316,29 +313,31 @@ Reflection *next_found_refl(Reflection *refl)
/********************************** Getters ***********************************/
+
/**
- * get_excitation_error:
+ * get_detector_pos:
* @refl: A %Reflection
+ * @fs: Location at which to store the fast scan offset of the reflection
+ * @ss: Location at which to store the slow scan offset of the reflection
*
- * Returns: The excitation error for the reflection.
**/
-double get_excitation_error(const Reflection *refl)
+void get_detector_pos(const Reflection *refl, double *fs, double *ss)
{
- return refl->data.excitation_error;
+ *fs = refl->data.fs;
+ *ss = refl->data.ss;
}
/**
- * get_detector_pos:
+ * get_panel:
* @refl: A %Reflection
- * @fs: Location at which to store the fast scan offset of the reflection
- * @ss: Location at which to store the slow scan offset of the reflection
+ *
+ * Returns: the panel which the reflection appears on
*
**/
-void get_detector_pos(const Reflection *refl, double *fs, double *ss)
+struct panel *get_panel(const Reflection *refl)
{
- *fs = refl->data.fs;
- *ss = refl->data.ss;
+ return refl->data.panel;
}
@@ -570,20 +569,32 @@ void copy_data(Reflection *to, const Reflection *from)
/**
* set_detector_pos:
* @refl: A %Reflection
- * @exerr: The excitation error for this reflection
* @fs: The fast scan offset of the reflection
* @ss: The slow scan offset of the reflection
*
**/
-void set_detector_pos(Reflection *refl, double exerr, double fs, double ss)
+void set_detector_pos(Reflection *refl, double fs, double ss)
{
- refl->data.excitation_error = exerr;
refl->data.fs = fs;
refl->data.ss = ss;
}
/**
+ * set_panel:
+ * @refl: A %Reflection
+ * @panel: Pointer to the panel structure on which the reflection appears
+ *
+ * Note that the pointer will be stored, not the contents of the structure.
+ *
+ **/
+void set_panel(Reflection *refl, struct panel *p)
+{
+ refl->data.panel = p;
+}
+
+
+/**
* set_partial:
* @refl: A %Reflection
* @rlow: The "low" excitation error
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index 85a87c54..8b0d8e04 100644
--- a/libcrystfel/src/reflist.h
+++ b/libcrystfel/src/reflist.h
@@ -3,11 +3,11 @@
*
* Fast reflection/peak list
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2011-2014 Thomas White <taw@physics.org>
+ * 2011-2015 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -84,8 +84,8 @@ extern Reflection *find_refl(const RefList *list, signed int h, signed int k, si
extern Reflection *next_found_refl(Reflection *refl);
/* Get */
-extern double get_excitation_error(const Reflection *refl);
extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
+extern struct panel *get_panel(const Reflection *refl);
extern double get_partiality(const Reflection *refl);
extern double get_lorentz(const Reflection *refl);
extern void get_indices(const Reflection *refl,
@@ -106,8 +106,8 @@ extern double get_mean_bg(const Reflection *refl);
/* Set */
extern void copy_data(Reflection *to, const Reflection *from);
-extern void set_detector_pos(Reflection *refl, double exerr,
- double fs, double ss);
+extern void set_detector_pos(Reflection *refl, double fs, double ss);
+extern void set_panel(Reflection *refl, struct panel *p);
extern void set_partial(Reflection *refl, double rlow, double rhigh, double p);
extern void set_partiality(Reflection *refl, double p);
extern void set_lorentz(Reflection *refl, double L);
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 8586c850..78ac9a40 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -319,7 +319,7 @@ static RefList *read_stream_reflections_2_3(FILE *fh, struct detector *det)
p = find_panel_by_name(det,pn);
write_fs = fs - p->orig_min_fs + p->min_fs;
write_ss = ss - p->orig_min_ss + p->min_ss;
- set_detector_pos(refl, 0.0, write_fs, write_ss);
+ set_detector_pos(refl, write_fs, write_ss);
}
set_esd_intensity(refl, sigma);
set_peak(refl, pk);
@@ -383,11 +383,11 @@ static RefList *read_stream_reflections_2_1(FILE *fh, struct detector *det)
p = find_orig_panel(det, fs, ss);
write_fs = fs - p->orig_min_fs + p->min_fs;
write_ss = ss - p->orig_min_ss + p->min_ss;
- set_detector_pos(refl, 0.0, write_fs, write_ss);
+ set_detector_pos(refl, write_fs, write_ss);
} else {
- set_detector_pos(refl, 0.0, fs, ss);
+ set_detector_pos(refl, fs, ss);
}
set_esd_intensity(refl, sigma);
@@ -448,11 +448,11 @@ static RefList *read_stream_reflections_2_2(FILE *fh, struct detector *det)
p = find_orig_panel(det, fs, ss);
write_fs = fs - p->orig_min_fs + p->min_fs;
write_ss = ss - p->orig_min_ss + p->min_ss;
- set_detector_pos(refl, 0.0, write_fs, write_ss);
+ set_detector_pos(refl, write_fs, write_ss);
} else {
- set_detector_pos(refl, 0.0, fs, ss);
+ set_detector_pos(refl, fs, ss);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-28 09:58:25 +0000