aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-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
-rw-r--r--src/partialator.c7
-rw-r--r--src/post-refinement.c20
-rw-r--r--src/predict-refine.c28
-rw-r--r--tests/integration_check.c2
-rw-r--r--tests/prediction_gradient_check.c20
12 files changed, 156 insertions, 256 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);
}
diff --git a/src/partialator.c b/src/partialator.c
index 4b1d5b18..ca3a9976 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -626,17 +626,12 @@ int main(int argc, char *argv[])
for ( j=0; j<images[i].n_crystals; j++ ) {
Crystal *cryst;
- int n_gained = 0;
- int n_lost = 0;
- double mean_p_change = 0.0;
cryst = images[i].crystals[j];
crystal_set_image(cryst, &images[i]);
/* Now it's safe to do the following */
- update_partialities_2(cryst, pmodel, &n_gained, &n_lost,
- &mean_p_change);
- assert(n_gained == 0); /* That'd just be silly */
+ update_partialities(cryst, pmodel);
}
}
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 96c0de12..745e5bd7 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -614,33 +614,17 @@ struct prdata pr_refine(Crystal *cr, const RefList *full,
double bsx, bsy, bsz;
double csx, csy, csz;
double dev;
- int n_total;
- int n_gained = 0;
- int n_lost = 0;
- n_total = num_reflections(crystal_get_reflections(cr));
cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
&bsx, &bsy, &bsz, &csx, &csy, &csz);
pr_iterate(cr, full, pmodel, &prdata.n_filtered);
- update_partialities_2(cr, pmodel, &n_gained, &n_lost,
- &mean_p_change);
+ update_partialities(cr, pmodel);
if ( verbose ) {
dev = guide_dev(cr, full);
- STATUS("PR Iteration %2i: mean p change = %10.2f"
- " dev = %10.5e, %i gained, %i lost, %i total\n",
- i+1, mean_p_change, dev,
- n_gained, n_lost, n_total);
- }
-
- if ( 3*n_lost > n_total ) {
- revert_crystal(cr, backup);
- update_partialities_2(cr, pmodel, &n_gained, &n_lost,
- &mean_p_change);
- crystal_set_user_flag(cr, 4);
- break;
+ STATUS("PR Iteration %2i: dev = %10.5e\n", i+1, dev);
}
i++;
diff --git a/src/predict-refine.c b/src/predict-refine.c
index 0e00dc5c..56a4cc2c 100644
--- a/src/predict-refine.c
+++ b/src/predict-refine.c
@@ -51,10 +51,12 @@ struct reflpeak {
Reflection *refl;
struct imagefeature *peak;
double Ih; /* normalised */
+ struct panel *panel; /* panel the reflection appears on
+ * (we assume this never changes) */
};
static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist,
- struct reflpeak *rps)
+ struct reflpeak *rps, struct detector *det)
{
int i;
const double min_dist = 0.25;
@@ -102,6 +104,7 @@ static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist,
rps[n_acc].refl = refl;
rps[n_acc].peak = f;
+ rps[n_acc].panel = find_panel(det, f->fs, f->ss);
n_acc++;
}
@@ -112,10 +115,9 @@ static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist,
static void twod_mapping(double fs, double ss, double *px, double *py,
- struct detector *det)
+ struct panel *p)
{
double xs, ys;
- struct panel *p = find_panel(det, fs, ss);
xs = fs*p->fsx + ss*p->ssx;
ys = fs*p->fsy + ss*p->ssy;
@@ -216,13 +218,13 @@ static double pos_gradient(int param, struct reflpeak *rp, struct detector *det,
double fsh, ssh;
double tt, clen, azi, azf;
- twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, det);
+ twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel);
get_detector_pos(rp->refl, &fsh, &ssh);
- twod_mapping(fsh, ssh, &xh, &yh, det);
+ twod_mapping(fsh, ssh, &xh, &yh, rp->panel);
get_indices(rp->refl, &h, &k, &l);
tt = asin(lambda * resolution(cell, h, k, l));
- clen = find_panel(det, fsh, ssh)->clen;
+ clen = rp->panel->clen;
azi = atan2(yh, xh);
azf = 2.0*(cos(azi) + sin(azi)); /* FIXME: Why factor of 2? */
@@ -276,9 +278,9 @@ static double pos_dev(struct reflpeak *rp, struct detector *det)
/* Peak position term */
double xpk, ypk, xh, yh;
double fsh, ssh;
- twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, det);
+ twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel);
get_detector_pos(rp->refl, &fsh, &ssh);
- twod_mapping(fsh, ssh, &xh, &yh, det);
+ twod_mapping(fsh, ssh, &xh, &yh, rp->panel);
return (xh-xpk) + (yh-ypk);
}
@@ -427,7 +429,7 @@ int refine_prediction(struct image *image, Crystal *cr)
if ( rps == NULL ) return 1;
n = pair_peaks(image->features, crystal_get_cell(cr),
- crystal_get_reflections(cr), rps);
+ crystal_get_reflections(cr), rps, image->det);
STATUS("%i peaks\n", n);
if ( n < 10 ) {
ERROR("Too few paired peaks to refine orientation.\n");
@@ -453,16 +455,10 @@ int refine_prediction(struct image *image, Crystal *cr)
/* Refine */
STATUS("Initial residual = %e\n", residual(rps, n, image->det));
for ( i=0; i<MAX_CYCLES; i++ ) {
- int n_gain = 0;
- int n_lost = 0;
- double mpc;
iterate(rps, n, crystal_get_cell(cr), image);
- update_partialities_2(cr, PMODEL_SCSPHERE, &n_gain, &n_lost,
- &mpc);
+ update_partialities(cr, PMODEL_SCSPHERE);
STATUS("Residual after iteration %i = %e\n",
i, residual(rps, n, image->det));
- STATUS("%i gained, %i lost, mean p change = %e\n", n_gain,
- n_lost, mpc);
}
free(rps);
diff --git a/tests/integration_check.c b/tests/integration_check.c
index 2d0a6c7a..ccb613c3 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -120,7 +120,7 @@ int main(int argc, char *argv[])
list = reflist_new();
refl = add_refl(list, 0, 0, 0);
- set_detector_pos(refl, 0.0, 64, 64);
+ set_detector_pos(refl, 64, 64);
cell = cell_new();
cell_set_lattice_type(cell, L_CUBIC);
cell_set_centering(cell, 'P');
diff --git a/tests/prediction_gradient_check.c b/tests/prediction_gradient_check.c
index 09b4a868..15b4d803 100644
--- a/tests/prediction_gradient_check.c
+++ b/tests/prediction_gradient_check.c
@@ -48,8 +48,7 @@ int checkrxy;
static void scan(RefList *reflections, RefList *compare,
- int *valid, long double *vals[3], int idx,
- struct detector *det)
+ int *valid, long double *vals[3], int idx)
{
int i;
Reflection *refl;
@@ -64,6 +63,7 @@ static void scan(RefList *reflections, RefList *compare,
Reflection *refl2;
double rlow, rhigh, p;
double fs, ss, xh, yh;
+ struct panel *panel;
get_indices(refl, &h, &k, &l);
refl2 = find_refl(compare, h, k, l);
@@ -75,7 +75,8 @@ static void scan(RefList *reflections, RefList *compare,
get_partial(refl2, &rlow, &rhigh, &p);
get_detector_pos(refl2, &fs, &ss);
- twod_mapping(fs, ss, &xh, &yh, det);
+ panel = get_panel(refl2);
+ twod_mapping(fs, ss, &xh, &yh, panel);
switch ( checkrxy ) {
@@ -169,16 +170,14 @@ static void calc_either_side(Crystal *cr, double incr_val,
cr_new = new_shifted_crystal(cr, refine, -incr_val);
compare = find_intersections(image, cr_new, PMODEL_SCSPHERE);
- scan(crystal_get_reflections(cr), compare, valid, vals, 0,
- crystal_get_image(cr)->det);
+ scan(crystal_get_reflections(cr), compare, valid, vals, 0);
cell_free(crystal_get_cell(cr_new));
crystal_free(cr_new);
reflist_free(compare);
cr_new = new_shifted_crystal(cr, refine, +incr_val);
compare = find_intersections(image, cr_new, PMODEL_SCSPHERE);
- scan(crystal_get_reflections(cr), compare, valid, vals, 2,
- crystal_get_image(cr)->det);
+ scan(crystal_get_reflections(cr), compare, valid, vals, 2);
cell_free(crystal_get_cell(cr_new));
crystal_free(cr_new);
reflist_free(compare);
@@ -231,8 +230,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
}
for ( i=0; i<nref; i++ ) valid[i] = 1;
- scan(reflections, reflections, valid, vals, 1,
- crystal_get_image(cr)->det);
+ scan(reflections, reflections, valid, vals, 1);
calc_either_side(cr, incr_val, valid, vals, refine);
@@ -282,11 +280,15 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
} else {
struct imagefeature pk;
+ struct image *image;
pk.fs = 0.0;
pk.ss = 0.0;
+
+ image = crystal_get_image(cr);
rp.refl = refl;
rp.peak = &pk;
+ rp.panel = &image->det->panels[0];
cgrad = pos_gradient(refine, &rp,
crystal_get_image(cr)->det,
generated by cgit v1.2.3 (git 2.25.1) at 2024-12-26 07:17:49 +0000