aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/calibrate_detector.c3
-rw-r--r--src/cubeit.c3
-rw-r--r--src/geometry.c52
-rw-r--r--src/geometry.h8
-rw-r--r--src/image.h33
-rw-r--r--src/indexamajig.c8
-rw-r--r--src/partialator.c81
-rw-r--r--src/pattern_sim.c2
-rw-r--r--src/peaks.c140
-rw-r--r--src/post-refinement.c127
-rw-r--r--src/post-refinement.h5
-rw-r--r--src/reflist.c139
-rw-r--r--src/reflist.h69
-rw-r--r--src/reintegrate.c3
-rw-r--r--src/templates.c52
15 files changed, 419 insertions, 306 deletions
diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c
index 2252842f..9aeae803 100644
--- a/src/calibrate_detector.c
+++ b/src/calibrate_detector.c
@@ -168,8 +168,7 @@ static void add_image(void *args, int cookie)
image.flags = NULL;
image.indexed_cell = NULL;
image.filename = pargs->filename;
- image.cpeaks = NULL;
- image.n_cpeaks = 0;
+ image.reflections = NULL;
image.det = NULL;
STATUS("%3i: Processing '%s'\n", cookie, pargs->filename);
diff --git a/src/cubeit.c b/src/cubeit.c
index 808e2e23..7b40f8e5 100644
--- a/src/cubeit.c
+++ b/src/cubeit.c
@@ -218,8 +218,7 @@ static void sum_image(void *pg, int cookie)
image.flags = NULL;
image.indexed_cell = NULL;
image.filename = apargs->filename;
- image.cpeaks = NULL;
- image.n_cpeaks = 0;
+ image.reflections = NULL;
image.det = pargs->det;
STATUS("Processing '%s'\n", apargs->filename);
diff --git a/src/geometry.c b/src/geometry.c
index 9bb41d12..9f9b953c 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -22,6 +22,7 @@
#include "image.h"
#include "peaks.h"
#include "beam-parameters.h"
+#include "reflist.h"
#define MAX_CPEAKS (256 * 256)
@@ -117,7 +118,7 @@ static double partiality(double r1, double r2, double r)
static int check_reflection(struct image *image, double mres, int output,
- struct cpeak *cpeaks, int np,
+ RefList *reflections,
signed int h, signed int k, signed int l,
double asx, double asy, double asz,
double bsx, double bsy, double bsz,
@@ -136,6 +137,7 @@ static int check_reflection(struct image *image, double mres, int output,
double divergence = image->div;
double lambda = image->lambda;
double klow, kcen, khigh; /* Wavenumber */
+ Reflection *refl;
/* "low" gives the largest Ewald sphere,
* "high" gives the smallest Ewald sphere. */
@@ -207,17 +209,9 @@ static int check_reflection(struct image *image, double mres, int output,
if ( p == -1 ) return 0;
/* Add peak to list */
- cpeaks[np].h = h;
- cpeaks[np].k = k;
- cpeaks[np].l = l;
- cpeaks[np].x = xda;
- cpeaks[np].y = yda;
- cpeaks[np].r1 = rlow;
- cpeaks[np].r2 = rhigh;
- cpeaks[np].p = part;
- cpeaks[np].clamp1 = clamp_low;
- cpeaks[np].clamp2 = clamp_high;
- np++;
+ refl = add_refl(reflections, h, k, l);
+ set_detector_pos(refl, 0.0, xda, yda);
+ set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
if ( output ) {
printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n",
@@ -228,23 +222,18 @@ static int check_reflection(struct image *image, double mres, int output,
}
-struct cpeak *find_intersections(struct image *image, UnitCell *cell,
- int *n, int output)
+RefList *find_intersections(struct image *image, UnitCell *cell,
+ int output)
{
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
- struct cpeak *cpeaks;
- int np = 0;
+ RefList *reflections;
int hmax, kmax, lmax;
double mres;
signed int h, k, l;
- cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
- if ( cpeaks == NULL ) {
- *n = 0;
- return NULL;
- }
+ reflections = reflist_new();
cell_get_reciprocal(cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
@@ -260,29 +249,30 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
for ( l=-lmax; l<lmax; l++ ) {
/* Ignore central beam */
if ( (h==0) && (k==0) && (l==0) ) continue;
- np += check_reflection(image, mres, output, cpeaks, np, h, k, l,
- asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
- if ( np == MAX_CPEAKS ) goto out;
+ check_reflection(image, mres, output, reflections, h, k, l,
+ asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
}
}
}
-out:
- *n = np;
- return cpeaks;
+ return reflections;
}
-double integrate_all(struct image *image, struct cpeak *cpeaks, int n)
+double integrate_all(struct image *image, RefList *reflections)
{
double itot = 0.0;
- int i;
+ Reflection *refl;
- for ( i=0; i<n; i++ ) {
+ for ( refl = first_refl(reflections);
+ refl != NULL;
+ refl = next_refl(refl) ) {
float x, y, intensity;
+ double xp, yp;
+ get_detector_pos(refl, &xp, &yp);
- if ( integrate_peak(image, cpeaks[i].x, cpeaks[i].y, &x, &y,
+ if ( integrate_peak(image, xp, yp, &x, &y,
&intensity, NULL, NULL, 0, 0) ) continue;
itot += intensity;
diff --git a/src/geometry.h b/src/geometry.h
index 0782a3e1..a0761456 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -17,10 +17,12 @@
#include <config.h>
#endif
-extern struct cpeak *find_intersections(struct image *image, UnitCell *cell,
- int *n, int output);
+#include "reflist.h"
-extern double integrate_all(struct image *image, struct cpeak *cpeaks, int n);
+extern RefList *find_intersections(struct image *image, UnitCell *cell,
+ int output);
+
+extern double integrate_all(struct image *image, RefList *reflections);
#endif /* GEOMETRY_H */
diff --git a/src/image.h b/src/image.h
index 10fb2406..ce0e7343 100644
--- a/src/image.h
+++ b/src/image.h
@@ -24,6 +24,7 @@
#include "utils.h"
#include "cell.h"
#include "detector.h"
+#include "reflist.h"
#define MAX_CELL_CANDIDATES (32)
@@ -52,35 +53,6 @@ struct imagefeature {
typedef struct _imagefeaturelist ImageFeatureList;
-/* This structure represents a predicted peak in an image */
-struct cpeak
-{
- /* Indices */
- signed int h;
- signed int k;
- signed int l;
-
- double min_distance;
- int valid;
-
- /* Partiality */
- double r1; /* First excitation error */
- double r2; /* Second excitation error */
- double p; /* Partiality */
- int clamp1; /* Clamp status for r1 */
- int clamp2; /* Clamp status for r2 */
-
- /* Location in image */
- int x;
- int y;
-
- int scalable;
-
- /* Intensity */
- double intensity;
-};
-
-
/* Structure describing an image */
struct image {
@@ -93,8 +65,7 @@ struct image {
struct detector *det;
struct beam_params *beam; /* The nominal beam parameters */
char *filename;
- struct cpeak *cpeaks;
- int n_cpeaks;
+ RefList *reflections;
int id; /* ID number of the thread
* handling this image */
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 6c01f95e..40d7cbd8 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -263,8 +263,7 @@ static struct image *get_simage(struct image *template, int alternate)
image->f0 = template->f0;
/* Prevent muppetry */
- image->cpeaks = NULL;
- image->n_cpeaks = 0;
+ image->reflections = NULL;
return image;
}
@@ -326,8 +325,7 @@ static void process_image(void *pp, int cookie)
image.indexed_cell = NULL;
image.id = cookie;
image.filename = filename;
- image.cpeaks = NULL;
- image.n_cpeaks = 0;
+ image.reflections = NULL;
image.det = pargs->static_args.det;
STATUS("Processing '%s'\n", image.filename);
@@ -434,7 +432,7 @@ done:
free(image.data);
free(image.flags);
image_feature_list_free(image.features);
- free(image.cpeaks);
+ reflist_free(image.reflections);
hdfile_close(hdfile);
}
diff --git a/src/partialator.c b/src/partialator.c
index d1ad3768..a4f3381c 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -34,6 +34,7 @@
#include "beam-parameters.h"
#include "post-refinement.h"
#include "hrs-scaling.h"
+#include "reflist.h"
/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
@@ -84,9 +85,9 @@ static void refine_image(int mytask, void *tasks)
struct image *image = pargs->image;
double nominal_photon_energy = pargs->image->beam->photon_energy;
struct hdfile *hdfile;
- struct cpeak *spots;
- int n, i;
+ int i;
double dev, last_dev;
+ RefList *reflections;
hdfile = hdfile_open(image->filename);
if ( hdfile == NULL ) {
@@ -110,16 +111,17 @@ static void refine_image(int mytask, void *tasks)
a/1.0e-9, b/1.0e-9, c/1.0e-9,
rad2deg(al), rad2deg(be), rad2deg(ga));
- spots = find_intersections(image, image->indexed_cell, &n, 0);
+ /* FIXME: Don't do this each time */
+ reflections = find_intersections(image, image->indexed_cell, 0);
dev = +INFINITY;
i = 0;
do {
last_dev = dev;
- dev = pr_iterate(image, pargs->i_full, pargs->sym, &spots, &n);
+ dev = pr_iterate(image, pargs->i_full, pargs->sym, reflections);
STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev);
i++;
} while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) );
- mean_partial_dev(image, spots, n, pargs->sym,
+ mean_partial_dev(image, reflections, pargs->sym,
pargs->i_full, pargs->graph);
if ( pargs->pgraph ) {
fprintf(pargs->pgraph, "%5i %5.2f\n", mytask, dev);
@@ -128,7 +130,7 @@ static void refine_image(int mytask, void *tasks)
free(image->data);
if ( image->flags != NULL ) free(image->flags);
hdfile_close(hdfile);
- free(spots);
+ reflist_free(reflections);
/* Muppet proofing */
image->data = NULL;
@@ -163,22 +165,23 @@ static void refine_all(struct image *images, int n_total_patterns,
}
-static void uniquify(struct cpeak *spot, const char *sym)
+static void uniquify(Reflection *refl, const char *sym)
{
+ signed int h, k, l;
signed int ha, ka, la;
- get_asymm(spot->h, spot->k, spot->l, &ha, &ka, &la, sym);
- spot->h = ha;
- spot->k = ka;
- spot->l = la;
+ get_indices(refl, &h, &k, &l);
+ get_asymm(h, k, l, &ha, &ka, &la, sym);
+ set_indices(refl, h, k, l);
}
+/* FIXME: Get rid of this */
static void integrate_image(struct image *image, ReflItemList *obs,
const char *sym)
{
- struct cpeak *spots;
- int j, n;
+ RefList *reflections;
+ Reflection *refl;
struct hdfile *hdfile;
double nominal_photon_energy = image->beam->photon_energy;
@@ -199,42 +202,39 @@ static void integrate_image(struct image *image, ReflItemList *obs,
}
/* Figure out which spots should appear in this pattern */
- spots = find_intersections(image, image->indexed_cell, &n, 0);
+ reflections = find_intersections(image, image->indexed_cell, 0);
/* For each reflection, estimate the partiality */
- for ( j=0; j<n; j++ ) {
+ for ( refl = first_refl(reflections);
+ refl != NULL;
+ refl = next_refl(refl) ) {
signed int h, k, l;
float i_partial;
float xc, yc;
+ double x, y;
- uniquify(&spots[j], sym);
-
- h = spots[j].h;
- k = spots[j].k;
- l = spots[j].l;
+ uniquify(refl, sym);
+ get_indices(refl, &h, &k, &l);
/* Don't attempt to use spots with very small
* partialities, since it won't be accurate. */
- if ( spots[j].p < 0.1 ) continue;
+ if ( get_partiality(refl) < 0.1 ) continue;
- /* Actual measurement of this reflection from this
- * pattern? */
- /* FIXME: Coordinates aren't whole numbers */
- if ( integrate_peak(image, spots[j].x, spots[j].y,
+ /* Actual measurement of this reflection from this pattern? */
+ get_detector_pos(refl, &x, &y);
+ if ( integrate_peak(image, x, y,
&xc, &yc, &i_partial, NULL, NULL, 1, 0) ) {
- spots[j].valid = 0;
+ delete_refl(refl);
continue;
}
- spots[j].intensity = i_partial;
- spots[j].valid = 1;
+ set_int(refl, i_partial);
if ( !find_item(obs, h, k, l) ) add_item(obs, h, k, l);
}
- image->cpeaks = spots;
- image->n_cpeaks = n;
+ image->reflections = reflections;
free(image->data);
if ( image->flags != NULL ) free(image->flags);
@@ -253,21 +253,20 @@ static void select_scalable_reflections(struct image *images, int n)
for ( m=0; m<n; m++ ) {
- int j;
+ Reflection *refl;
- for ( j=0; j<images[m].n_cpeaks; j++ ) {
+ for ( refl = first_refl(images[m].reflections);
+ refl != NULL;
+ refl = next_refl(refl) ) {
int scalable = 1;
+ double v;
- if ( images[m].cpeaks[j].p < 0.1 ) scalable = 0;
- if ( !images[m].cpeaks[j].valid ) {
- scalable = 0;
- } else {
- double v = fabs(images[m].cpeaks[j].intensity);
- if ( v < 0.1 ) scalable = 0;
- }
+ if ( get_partiality(refl) < 0.1 ) scalable = 0;
+ v = fabs(get_intensity(refl));
+ if ( v < 0.1 ) scalable = 0;
- images[m].cpeaks[j].scalable = scalable;
+ set_scalable(refl, scalable);
}
@@ -525,7 +524,7 @@ int main(int argc, char *argv[])
/* Clean up */
for ( i=0; i<n_total_patterns; i++ ) {
- free(images[i].cpeaks);
+ reflist_free(images[i].reflections);
}
free(I_full);
delete_items(obs);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 4584bb3c..90425328 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -567,7 +567,7 @@ int main(int argc, char *argv[])
find_projected_peaks(&image, cell, 0, 0.1);
output_intensities(&image, cell, NULL, 0, 0, stdout,
0, 0.1);
- free(image.cpeaks);
+ reflist_free(image.reflections);
}
if ( powder_fn != NULL ) {
diff --git a/src/peaks.c b/src/peaks.c
index 88ac50bf..d1132a3a 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -32,9 +32,6 @@
#include "diffraction.h"
-/* Maximum number of peaks which may be predicted by find_projected_peaks() */
-#define MAX_CPEAKS (8192)
-
/* How close a peak must be to an indexed position to be considered "close"
* for the purposes of double hit detection and sanity checking. */
#define PEAK_CLOSE (30.0)
@@ -483,12 +480,11 @@ int find_projected_peaks(struct image *image, UnitCell *cell,
double ax, ay, az;
double bx, by, bz;
double cx, cy, cz;
- struct cpeak *cpeaks;
- int n_cpeaks = 0;
+ RefList *reflections;
double alen, blen, clen;
+ int n_reflections = 0;
- cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
- if ( cpeaks == NULL ) return 0;
+ reflections = reflist_new();
/* "Borrow" direction values to get reciprocal lengths */
cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
@@ -507,8 +503,8 @@ int find_projected_peaks(struct image *image, UnitCell *cell,
signed int h, k, l;
struct rvec q;
double dist;
- int found = 0;
- int j;
+ Reflection *refl;
+ double cur_dist;
q = get_q(image, x, y, 1, NULL, 1.0/image->lambda);
@@ -535,44 +531,26 @@ int find_projected_peaks(struct image *image, UnitCell *cell,
if ( dist > domain_r ) continue;
}
- for ( j=0; j<n_cpeaks; j++ ) {
- if ( (cpeaks[j].h == h) && (cpeaks[j].k == k)
- && (cpeaks[j].l == l) ) {
-
- if ( dist < cpeaks[j].min_distance ) {
- cpeaks[j].min_distance = dist;
- cpeaks[j].x = x;
- cpeaks[j].y = y;
- }
-
- found = 1;
-
- }
- }
-
- if ( !found ) {
- cpeaks[n_cpeaks].min_distance = dist;
- cpeaks[n_cpeaks].x = x;
- cpeaks[n_cpeaks].y = y;
- cpeaks[n_cpeaks].h = h;
- cpeaks[n_cpeaks].k = k;
- cpeaks[n_cpeaks].l = l;
- n_cpeaks++;
- if ( n_cpeaks == MAX_CPEAKS ) {
- ERROR("Unit cell is insanely large!\n");
- goto out;
+ refl = find_refl(reflections, h, k, l);
+ if ( refl != NULL ) {
+ cur_dist = get_excitation_error(refl);
+ if ( dist < cur_dist ) {
+ set_detector_pos(refl, dist, x, y);
}
+ } else {
+ add_refl_with_det_pos(reflections, h, k, l, x, y, dist);
+ n_reflections++;
}
}
}
-out:
- STATUS("Found %i reflections\n", n_cpeaks);
- image->cpeaks = cpeaks;
- image->n_cpeaks = n_cpeaks;
+ optimise_reflist(reflections);
+
+ STATUS("Found %i reflections\n", n_reflections);
+ image->reflections = reflections;
- return n_cpeaks;
+ return n_reflections;
}
@@ -685,20 +663,14 @@ void output_intensities(struct image *image, UnitCell *cell,
int use_closer, FILE *ofh,
int circular_domain, double domain_r)
{
- int i;
- int n_found;
- int n_indclose = 0;
- int n_foundclose = 0;
- int n_veto = 0;
- int n_veto_second = 0;
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
+ Reflection *refl;
- if ( image->n_cpeaks == 0 ) {
+ if ( image->reflections == NULL ) {
find_projected_peaks(image, cell, circular_domain, domain_r);
}
- if ( image->n_cpeaks == 0 ) return;
/* Get exclusive access to the output stream if necessary */
if ( mutex != NULL ) pthread_mutex_lock(mutex);
@@ -709,16 +681,20 @@ void output_intensities(struct image *image, UnitCell *cell,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
- for ( i=0; i<image->n_cpeaks; i++ ) {
+ for ( refl = first_refl(image->reflections);
+ refl != NULL;
+ refl = next_refl(refl) ) {
float x, y, intensity;
double d;
int idx;
double bg, max;
struct panel *p;
+ double px, py;
+ signed int h, k, l;
- p = find_panel(image->det, image->cpeaks[i].x,
- image->cpeaks[i].y);
+ get_detector_pos(refl, &px, &py);
+ p = find_panel(image->det, px, py);
if ( p == NULL ) continue;
if ( p->no_index ) continue;
@@ -729,9 +705,7 @@ void output_intensities(struct image *image, UnitCell *cell,
if ( image->features != NULL ) {
f = image_feature_closest(image->features,
- image->cpeaks[i].x,
- image->cpeaks[i].y,
- &d, &idx);
+ px, py, &d, &idx);
} else {
f = NULL;
}
@@ -754,7 +728,6 @@ void output_intensities(struct image *image, UnitCell *cell,
* region. Integration of the same
* peak included a bad region this time.
*/
- n_veto_second++;
continue;
}
intensity = f->intensity;
@@ -763,81 +736,36 @@ void output_intensities(struct image *image, UnitCell *cell,
int r;
- r = integrate_peak(image,
- image->cpeaks[i].x,
- image->cpeaks[i].y,
- &x, &y, &intensity, &bg, &max,
+ r = integrate_peak(image, px, py, &x, &y,
+ &intensity, &bg, &max,
polar, 1);
if ( r ) {
/* Plain old ordinary peak veto */
- n_veto++;
continue;
}
}
- if ( (f != NULL) && (d < PEAK_CLOSE) ) {
- n_indclose++;
- }
-
} else {
int r;
- r = integrate_peak(image,
- image->cpeaks[i].x,
- image->cpeaks[i].y,
- &x, &y, &intensity, &bg, &max, polar,
- 0);
+ r = integrate_peak(image, px, py, &x, &y,
+ &intensity, &bg, &max, polar, 0);
if ( r ) {
/* Plain old ordinary peak veto */
- n_veto++;
continue;
}
}
/* Write h,k,l, integrated intensity and centroid coordinates */
+ get_indices(refl, &h, &k, &l);
fprintf(ofh, "%3i %3i %3i %6f (at %5.2f,%5.2f) max=%6f bg=%6f\n",
- image->cpeaks[i].h, image->cpeaks[i].k, image->cpeaks[i].l,
- intensity, x, y, max, bg);
-
- image->cpeaks[i].x = x;
- image->cpeaks[i].y = y;
-
- }
-
- n_found = image_feature_count(image->features);
- for ( i=0; i<n_found; i++ ) {
-
- struct imagefeature *f;
- int j;
-
- f = image_get_feature(image->features, i);
-
- if ( f == NULL ) continue;
-
- for ( j=0; j<image->n_cpeaks; j++ ) {
-
- double d;
-
- d = pow(image->cpeaks[j].x-f->x, 2.0)
- + pow(image->cpeaks[j].y-f->y, 2.0);
-
- if ( d < PEAK_CLOSE ) n_foundclose++;
-
- }
+ h, l, l, intensity, x, y, max, bg);
}
- fprintf(ofh, "Peak statistics: %i peaks found by the peak search out of "
- "%i were close to indexed positions. "
- "%i indexed positions out of %i were close to detected peaks.\n",
- n_foundclose, n_found, n_indclose, image->n_cpeaks);
- fprintf(ofh, "%i integrations using indexed locations were aborted because "
- "they hit one or more bad pixels.\n", n_veto);
- fprintf(ofh, "%i integrations using peak search locations were aborted "
- "because they hit one or more bad pixels.\n", n_veto_second);
/* Blank line at end */
fprintf(ofh, "\n");
diff --git a/src/post-refinement.c b/src/post-refinement.c
index b860bb42..25664e69 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -65,7 +65,7 @@ static double partiality_rgradient(double r, double profile_radius)
/* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(struct image *image, int k, struct cpeak spot, double r)
+double gradient(struct image *image, int k, Reflection *refl, double r)
{
double ds, tt, azi;
double nom, den;
@@ -74,24 +74,32 @@ double gradient(struct image *image, int k, struct cpeak spot, double r)
double bsx, bsy, bsz;
double csx, csy, csz;
double xl, yl, zl;
+ signed int hi, ki, li;
+ double r1, r2, p;
+ int clamp_low, clamp_high;
+
+ get_indices(refl, &hi, &ki, &li);
cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
- xl = spot.h*asx + spot.k*bsx + spot.l*csx;
- yl = spot.h*asy + spot.k*bsy + spot.l*csy;
- zl = spot.h*asz + spot.k*bsz + spot.l*csz;
+ xl = hi*asx + ki*bsx + li*csx;
+ yl = hi*asy + ki*bsy + li*csy;
+ zl = hi*asz + ki*bsz + li*csz;
+
- ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l);
+ ds = 2.0 * resolution(image->indexed_cell, hi, ki, li);
tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+k);
azi = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0);
+ get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high);
+
/* Calculate the gradient of partiality wrt excitation error. */
- if ( spot.clamp1 == 0 ) {
- g += partiality_gradient(spot.r1, r);
+ if ( clamp_low == 0 ) {
+ g += partiality_gradient(r1, r);
}
- if ( spot.clamp2 == 0 ) {
- g += partiality_gradient(spot.r2, r);
+ if ( clamp_high == 0 ) {
+ g += partiality_gradient(r2, r);
}
/* For many gradients, just multiply the above number by the gradient
@@ -99,7 +107,7 @@ double gradient(struct image *image, int k, struct cpeak spot, double r)
switch ( k ) {
case REF_SCALE :
- return -spot.p*pow(image->osf, -2.0);
+ return -p*pow(image->osf, -2.0);
case REF_DIV :
nom = sqrt(2.0) * ds * sin(image->div/2.0);
@@ -107,33 +115,33 @@ double gradient(struct image *image, int k, struct cpeak spot, double r)
return (nom/den) * g;
case REF_R :
- if ( spot.clamp1 == 0 ) {
- g += partiality_rgradient(spot.r1, r);
+ if ( clamp_low == 0 ) {
+ g += partiality_rgradient(r1, r);
}
- if ( spot.clamp2 == 0 ) {
- g += partiality_rgradient(spot.r2, r);
+ if ( clamp_high == 0 ) {
+ g += partiality_rgradient(r2, r);
}
return g;
/* Cell parameters and orientation */
case REF_ASX :
- return spot.h * pow(sin(tt), -1.0) * cos(azi) * g;
+ return hi * pow(sin(tt), -1.0) * cos(azi) * g;
case REF_BSX :
- return spot.k * pow(sin(tt), -1.0) * cos(azi) * g;
+ return ki * pow(sin(tt), -1.0) * cos(azi) * g;
case REF_CSX :
- return spot.l * pow(sin(tt), -1.0) * cos(azi) * g;
+ return li * pow(sin(tt), -1.0) * cos(azi) * g;
case REF_ASY :
- return spot.h * pow(sin(tt), -1.0) * sin(azi) * g;
+ return hi * pow(sin(tt), -1.0) * sin(azi) * g;
case REF_BSY :
- return spot.k * pow(sin(tt), -1.0) * sin(azi) * g;
+ return ki * pow(sin(tt), -1.0) * sin(azi) * g;
case REF_CSY :
- return spot.l * pow(sin(tt), -1.0) * sin(azi) * g;
+ return li * pow(sin(tt), -1.0) * sin(azi) * g;
case REF_ASZ :
- return spot.h * pow(cos(tt), -1.0) * g;
+ return hi * pow(cos(tt), -1.0) * g;
case REF_BSZ :
- return spot.k * pow(cos(tt), -1.0) * g;
+ return ki * pow(cos(tt), -1.0) * g;
case REF_CSZ :
- return spot.l * pow(cos(tt), -1.0) * g;
+ return li * pow(cos(tt), -1.0) * g;
}
@@ -218,31 +226,34 @@ void apply_shift(struct image *image, int k, double shift)
}
-double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
+double mean_partial_dev(struct image *image, RefList *reflections,
const char *sym, double *i_full, FILE *graph)
{
- int h, n_used;
+ int n_used;
double delta_I = 0.0;
+ Reflection *refl;
n_used = 0;
- for ( h=0; h<n; h++ ) {
+ for ( refl = first_refl(reflections);
+ refl != NULL;
+ refl = next_refl(refl) ) {
signed int hind, kind, lind;
signed int ha, ka, la;
double I_full;
float I_partial;
float xc, yc;
+ double x, y;
+ double p;
- hind = spots[h].h;
- kind = spots[h].k;
- lind = spots[h].l;
+ get_indices(refl, &hind, &kind, &lind);
- if ( !spots[h].scalable ) continue;
+ if ( !get_scalable(refl) ) continue;
- /* Actual measurement of this reflection from this
- * pattern? */
- /* FIXME: Coordinates aren't whole numbers */
- if ( integrate_peak(image, spots[h].x, spots[h].y,
+ /* Actual measurement of this reflection from this pattern? */
+ get_detector_pos(refl, &x, &y);
+ /* FIXME: Get rid of this */
+ if ( integrate_peak(image, x, y,
&xc, &yc, &I_partial, NULL, NULL,
1, 0) ) {
continue;
@@ -250,14 +261,15 @@ double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
I_full = lookup_intensity(i_full, ha, ka, la);
- delta_I += fabs(I_partial - (spots[h].p * I_full / image->osf));
+ p = get_partiality(refl);
+ delta_I += fabs(I_partial - (p * I_full / image->osf));
n_used++;
if ( graph != NULL ) {
fprintf(graph, "%3i %3i %3i %5.2f (at %5.2f,%5.2f)"
" %5.2f %5.2f\n",
- hind, kind, lind, I_partial/spots[h].p, xc, yc,
- spots[h].p, I_partial / I_full);
+ hind, kind, lind, I_partial/p, xc, yc,
+ p, I_partial / I_full);
}
}
@@ -268,19 +280,21 @@ double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
/* Perform one cycle of post refinement on 'image' against 'i_full' */
double pr_iterate(struct image *image, double *i_full, const char *sym,
- struct cpeak **pspots, int *n)
+ RefList *reflections)
{
gsl_matrix *M;
gsl_vector *v;
gsl_vector *shifts;
- int h, param;
- struct cpeak *spots = *pspots;
+ int param;
+ Reflection *refl;
M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
v = gsl_vector_calloc(NUM_PARAMS);
/* Construct the equations, one per reflection in this image */
- for ( h=0; h<*n; h++ ) {
+ for ( refl = first_refl(reflections);
+ refl != NULL;
+ refl = next_refl(refl) ) {
signed int hind, kind, lind;
signed int ha, ka, la;
@@ -288,17 +302,17 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
float I_partial;
float xc, yc;
int k;
+ double x, y;
+ double p;
- hind = spots[h].h;
- kind = spots[h].k;
- lind = spots[h].l;
+ get_indices(refl, &hind, &kind, &lind);
- if ( !spots[h].scalable ) continue;
+ if ( !get_scalable(refl) ) continue;
- /* Actual measurement of this reflection from this
- * pattern? */
- /* FIXME: Coordinates aren't whole numbers */
- if ( integrate_peak(image, spots[h].x, spots[h].y,
+ /* Actual measurement of this reflection from this pattern? */
+ get_detector_pos(refl, &x, &y);
+ /* FIXME: Get rid of this */
+ if ( integrate_peak(image, x, y,
&xc, &yc, &I_partial, NULL, NULL,
1, 0) ) {
continue;
@@ -306,7 +320,8 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
I_full = lookup_intensity(i_full, ha, ka, la);
- delta_I = I_partial - (spots[h].p * I_full / image->osf);
+ p = get_partiality(refl);
+ delta_I = I_partial - (p * I_full / image->osf);
for ( k=0; k<NUM_PARAMS; k++ ) {
@@ -319,9 +334,9 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
M_curr = gsl_matrix_get(M, g, k);
- M_c = gradient(image, g, spots[h],
+ M_c = gradient(image, g, refl,
image->profile_radius)
- * gradient(image, k, spots[h],
+ * gradient(image, k, refl,
image->profile_radius);
M_c *= pow(I_full, 2.0);
@@ -329,8 +344,7 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
}
- gr = gradient(image, k, spots[h],
- image->profile_radius);
+ gr = gradient(image, k, refl, image->profile_radius);
v_c = delta_I * I_full * gr;
gsl_vector_set(v, k, v_c);
@@ -350,8 +364,5 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
gsl_vector_free(v);
gsl_vector_free(shifts);
- free(spots);
- spots = find_intersections(image, image->indexed_cell, n, 0);
- *pspots = spots;
- return mean_partial_dev(image, spots, *n, sym, i_full, NULL);
+ return mean_partial_dev(image, reflections, sym, i_full, NULL);
}
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 983ecc28..b3c6c074 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -21,6 +21,7 @@
#include <stdio.h>
#include "image.h"
+#include "reflist.h"
/* Refineable parameters */
@@ -44,12 +45,12 @@ enum {
void apply_shift(struct image *image, int k, double shift);
-double mean_partial_dev(struct image *image, struct cpeak *spots, int n,
+double mean_partial_dev(struct image *image, RefList *reflections,
const char *sym, double *i_full, FILE *graph);
double pr_iterate(struct image *image, double *i_full, const char *sym,
- struct cpeak **pspots, int *n);
+ RefList *reflections);
#endif /* POST_REFINEMENT_H */
diff --git a/src/reflist.c b/src/reflist.c
new file mode 100644
index 00000000..990bc535
--- /dev/null
+++ b/src/reflist.c
@@ -0,0 +1,139 @@
+/*
+ * reflist.c
+ *
+ * Fast reflection/peak list
+ *
+ * (c) 2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#include <stdlib.h>
+
+#include "reflist.h"
+
+
+struct _reflection {
+
+ unsigned int serial; /* Serial number */
+
+ signed int h;
+ signed int k;
+ signed int l;
+
+ double excitation_error;
+
+ /* Partiality */
+ double r1; /* First excitation error */
+ double r2; /* Second excitation error */
+ double p; /* Partiality */
+ int clamp1; /* Clamp status for r1 */
+ int clamp2; /* Clamp status for r2 */
+
+ /* Location in image */
+ int x;
+ int y;
+
+ int scalable;
+
+ /* Intensity */
+ double intensity;
+
+
+};
+
+
+struct _reflist {
+
+ struct ref *head;
+
+};
+
+
+/**************************** Creation / deletion *****************************/
+
+/* Create a reflection list */
+RefList *reflist_new()
+{
+ RefList *new;
+
+ new = malloc(sizeof(struct _reflist));
+
+ return new;
+}
+
+
+void reflist_free(RefList *list)
+{
+}
+
+
+/********************************** Search ************************************/
+
+Reflection *find_refl(RefList *list, INDICES)
+{
+}
+
+
+/********************************** Getters ***********************************/
+
+double get_excitation_error(Reflection *refl)
+{
+}
+
+
+void get_detector_pos(Reflection *refl, double *x, double *y)
+{
+}
+
+
+void get_indices(Reflection *refl, signed int *h, signed int *k, signed int *l)
+{
+}
+
+
+/********************************** Setters ***********************************/
+
+void set_detector_pos(Reflection *refl, double exerr, double x, double y)
+{
+}
+
+
+void set_partial(Reflection *refl, double r1, double r2, double p,
+ double clamp_low, double clamp_high)
+{
+}
+
+
+/********************************* Insertion **********************************/
+
+Reflection *add_refl(RefList *list, INDICES)
+{
+}
+
+
+Reflection *add_refl_with_det_pos(RefList *refl, INDICES, double exerr,
+ double x, double y)
+{
+}
+
+
+/********************************* Iteration **********************************/
+
+Reflection *first_refl(RefList *list)
+{
+}
+
+
+Reflection *next_refl(Reflection *refl)
+{
+}
+
+
+/*********************************** Voodoo ***********************************/
+
+void optimise_reflist(RefList *list)
+{
+}
diff --git a/src/reflist.h b/src/reflist.h
new file mode 100644
index 00000000..0ba09a3b
--- /dev/null
+++ b/src/reflist.h
@@ -0,0 +1,69 @@
+/*
+ * reflist.h
+ *
+ * Fast reflection/peak list
+ *
+ * (c) 2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifndef REFLIST_H
+#define REFLIST_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+typedef struct _reflist RefList;
+typedef struct _reflection Reflection;
+
+#define INDICES signed int h, signed int k, signed int l
+
+/* Creation/deletion */
+extern RefList *reflist_new(void);
+extern void reflist_free(RefList *list);
+
+/* Search */
+extern Reflection *find_refl(RefList *list, INDICES);
+
+/* Get */
+extern double get_excitation_error(Reflection *refl);
+extern void get_detector_pos(Reflection *refl, double *x, double *y);
+extern void get_indices(Reflection *refl,
+ signed int *h, signed int *k, signed int *l);
+extern double get_partiality(Reflection *refl);
+extern double get_intensity(Reflection *refl);
+extern void get_partial(Reflection *refl, double *r1, double *r2, double *p,
+ int *clamp_low, int *clamp_high);
+extern int get_scalable(Reflection *refl);
+
+/* Set */
+extern void set_detector_pos(Reflection *refl, double exerr,
+ double x, double y);
+extern void set_partial(Reflection *refl, double r1, double r2, double p,
+ double clamp_low, double clamp_high);
+extern void set_indices(Reflection *refl,
+ signed int h, signed int k, signed int l);
+extern void set_int(Reflection *refl, double intensity);
+extern void set_scalable(Reflection *refl, int scalable);
+
+/* Insertion */
+extern Reflection *add_refl(RefList *list, INDICES);
+extern Reflection *add_refl_with_det_pos(RefList *list, INDICES, double exerr,
+ double x, double y);
+
+
+/* Deletion */
+extern void delete_refl(Reflection *refl);
+
+/* Iteration */
+extern Reflection *first_refl(RefList *list);
+extern Reflection *next_refl(Reflection *refl);
+
+/* Voodoo */
+extern void optimise_reflist(RefList *list);
+
+#endif /* REFLIST_H */
diff --git a/src/reintegrate.c b/src/reintegrate.c
index 9d882d4b..09634cd7 100644
--- a/src/reintegrate.c
+++ b/src/reintegrate.c
@@ -110,8 +110,7 @@ static void process_image(void *pg, int cookie)
image.flags = NULL;
image.indexed_cell = NULL;
image.filename = apargs->filename;
- image.cpeaks = NULL;
- image.n_cpeaks = 0;
+ image.reflections = NULL;
image.det = pargs->det;
STATUS("Processing '%s'\n", apargs->filename);
diff --git a/src/templates.c b/src/templates.c
index 4bd85a76..9abbcb17 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -21,6 +21,7 @@
#include "geometry.h"
#include "hdf5-file.h"
#include "peaks.h"
+#include "reflist.h"
#include <assert.h>
@@ -42,8 +43,7 @@ struct _indexingprivate_template
struct template {
double omega;
double phi;
- int n;
- struct cpeak *spots;
+ RefList *spots;
};
@@ -164,24 +164,22 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename,
for ( omega = 0.0; omega < omega_max-omega_step; omega += omega_step ) {
for ( phi = 0.0; phi < phi_max-phi_step; phi += phi_step ) {
- int n;
- struct cpeak *cpeaks;
+ RefList *reflections;
UnitCell *cell_rot;
assert(i < n_templates);
cell_rot = rotate_cell(cell, omega, phi, 0.0);
- cpeaks = find_intersections(&image, cell_rot, &n, 0);
- if ( cpeaks == NULL ) {
+ reflections = find_intersections(&image, cell_rot, 0);
+ if ( reflections == NULL ) {
ERROR("Template calculation failed.\n");
return NULL;
}
priv->templates[i].omega = omega;
priv->templates[i].phi = phi;
- priv->templates[i].n = n;
- priv->templates[i].spots = cpeaks;
+ priv->templates[i].spots = reflections;
i++;
free(cell_rot);
@@ -227,23 +225,28 @@ static int fast_integrate_peak(struct image *image, int xp, int yp)
}
-static double integrate_all_rot(struct image *image, struct cpeak *cpeaks,
- int n, double rot)
+static double integrate_all_rot(struct image *image, RefList *reflections,
+ double rot)
{
double itot = 0.0;
- int i;
double cosr, sinr;
int num_int = 0;
+ Reflection *refl;
cosr = cos(rot);
sinr = sin(rot);
- for ( i=0; i<n; i++ ) {
+ for ( refl = first_refl(reflections);
+ refl != NULL;
+ refl = next_refl(refl) ) {
float xp, yp;
+ double x, y;
+
+ get_detector_pos(refl, &x, &y);
- xp = cosr*cpeaks[i].x + sinr*cpeaks[i].y;
- yp = -sinr*cpeaks[i].x + cosr*cpeaks[i].y;
+ xp = cosr*x + sinr*y;
+ yp = -sinr*x + cosr*y;
itot += fast_integrate_peak(image, rint(xp), rint(yp));
num_int++;
@@ -256,26 +259,31 @@ static double integrate_all_rot(struct image *image, struct cpeak *cpeaks,
/* Return the mean of the distances between peaks in the image and peaks from
* the given template. */
-static double mean_distance(struct image *image, struct cpeak *cpeaks,
- int n, double rot)
+static double mean_distance(struct image *image, RefList *reflections,
+ double rot)
{
double dtot = 0.0;
- int i;
double cosr, sinr;
int num_dist = 0;
+ Reflection *refl;
cosr = cos(rot);
sinr = sin(rot);
/* For each template peak */
- for ( i=0; i<n; i++ ) {
+ for ( refl = first_refl(reflections);
+ refl != NULL;
+ refl = next_refl(refl) ) {
float xp, yp;
int j;
double min_dsq;
+ double x, y;
+
+ get_detector_pos(refl, &x, &y);
- xp = cosr*cpeaks[i].x + sinr*cpeaks[i].y;
- yp = -sinr*cpeaks[i].x + cosr*cpeaks[i].y;
+ xp = cosr*x + sinr*y;
+ yp = -sinr*x + cosr*y;
/* Compare to every real peak */
min_dsq = +INFINITY;
@@ -328,11 +336,11 @@ void match_templates(struct image *image, IndexingPrivate *ipriv)
if ( !peaks ) {
val = integrate_all_rot(image, priv->templates[i].spots,
- priv->templates[i].n, rot);
+ rot);
best = val > max;
} else {
val = mean_distance(image, priv->templates[i].spots,
- priv->templates[i].n, rot);
+ rot);
best = val < max;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-08 17:57:45 +0000