aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/partialator.c242
-rw-r--r--src/post-refinement.c105
-rw-r--r--src/post-refinement.h33
-rw-r--r--src/reflist.c27
-rw-r--r--src/reflist.h2
5 files changed, 119 insertions, 290 deletions
diff --git a/src/partialator.c b/src/partialator.c
index a4f3381c..37c8fbd6 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -37,10 +37,6 @@
#include "reflist.h"
-/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
-#define MAX_CYCLES (100)
-
-
static void show_help(const char *s)
{
printf("Syntax: %s [options]\n\n", s);
@@ -57,9 +53,6 @@ static void show_help(const char *s)
" initial values for parameters, and nominal\n"
" wavelengths if no per-shot value is found in \n"
" an HDF5 file.\n"
-" -x, --prefix=<p> Prefix filenames from input file with <p>.\n"
-" --basename Remove the directory parts of the filenames.\n"
-" --no-check-prefix Don't attempt to correct the --prefix.\n"
" -y, --symmetry=<sym> Merge according to symmetry <sym>.\n"
" -n, --iterations=<n> Run <n> cycles of scaling and post-refinement.\n"
"\n"
@@ -83,58 +76,8 @@ static void refine_image(int mytask, void *tasks)
struct refine_args *all_args = tasks;
struct refine_args *pargs = &all_args[mytask];
struct image *image = pargs->image;
- double nominal_photon_energy = pargs->image->beam->photon_energy;
- struct hdfile *hdfile;
- int i;
- double dev, last_dev;
- RefList *reflections;
-
- hdfile = hdfile_open(image->filename);
- if ( hdfile == NULL ) {
- ERROR("Couldn't open '%s'\n", image->filename);
- return;
- } else if ( hdfile_set_image(hdfile, "/data/data0") ) {
- ERROR("Couldn't select path\n");
- hdfile_close(hdfile);
- return;
- }
-
- if ( hdf5_read(hdfile, pargs->image, 0, nominal_photon_energy) ) {
- ERROR("Couldn't read '%s'\n", image->filename);
- hdfile_close(hdfile);
- return;
- }
-
- double a, b, c, al, be, ga;
- cell_get_parameters(image->indexed_cell, &a, &b, &c, &al, &be, &ga);
- STATUS("Initial cell: %5.2f %5.2f %5.2f nm %5.2f %5.2f %5.2f deg\n",
- a/1.0e-9, b/1.0e-9, c/1.0e-9,
- rad2deg(al), rad2deg(be), rad2deg(ga));
-
- /* 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, 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, reflections, pargs->sym,
- pargs->i_full, pargs->graph);
- if ( pargs->pgraph ) {
- fprintf(pargs->pgraph, "%5i %5.2f\n", mytask, dev);
- }
- free(image->data);
- if ( image->flags != NULL ) free(image->flags);
- hdfile_close(hdfile);
- reflist_free(reflections);
-
- /* Muppet proofing */
- image->data = NULL;
- image->flags = NULL;
+ pr_refine(image, pargs->i_full, pargs->sym);
}
@@ -165,87 +108,6 @@ static void refine_all(struct image *images, int n_total_patterns,
}
-static void uniquify(Reflection *refl, const char *sym)
-{
- signed int h, k, l;
- signed int ha, ka, 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)
-{
- RefList *reflections;
- Reflection *refl;
- struct hdfile *hdfile;
- double nominal_photon_energy = image->beam->photon_energy;
-
- hdfile = hdfile_open(image->filename);
- if ( hdfile == NULL ) {
- ERROR("Couldn't open '%s'\n", image->filename);
- return;
- } else if ( hdfile_set_image(hdfile, "/data/data0") ) {
- ERROR("Couldn't select path\n");
- hdfile_close(hdfile);
- return;
- }
-
- if ( hdf5_read(hdfile, image, 0, nominal_photon_energy) ) {
- ERROR("Couldn't read '%s'\n", image->filename);
- hdfile_close(hdfile);
- return;
- }
-
- /* Figure out which spots should appear in this pattern */
- reflections = find_intersections(image, image->indexed_cell, 0);
-
- /* For each reflection, estimate the partiality */
- 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(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 ( get_partiality(refl) < 0.1 ) continue;
-
- /* 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) ) {
- delete_refl(refl);
- continue;
- }
-
- set_int(refl, i_partial);
-
- if ( !find_item(obs, h, k, l) ) add_item(obs, h, k, l);
-
- }
- image->reflections = reflections;
-
- free(image->data);
- if ( image->flags != NULL ) free(image->flags);
- hdfile_close(hdfile);
-
- /* Muppet proofing */
- image->data = NULL;
- image->flags = NULL;
-}
-
-
/* Decide which reflections can be scaled */
static void select_scalable_reflections(struct image *images, int n)
{
@@ -280,12 +142,9 @@ int main(int argc, char *argv[])
char *infile = NULL;
char *outfile = NULL;
char *geomfile = NULL;
- char *prefix = NULL;
char *sym = NULL;
FILE *fh;
int nthreads = 1;
- int config_basename = 0;
- int config_checkprefix = 1;
struct detector *det;
unsigned int *cts;
ReflItemList *obs;
@@ -303,9 +162,6 @@ int main(int argc, char *argv[])
{"output", 1, NULL, 'o'},
{"geometry", 1, NULL, 'g'},
{"beam", 1, NULL, 'b'},
- {"prefix", 1, NULL, 'x'},
- {"basename", 0, &config_basename, 1},
- {"no-check-prefix", 0, &config_checkprefix, 0},
{"symmetry", 1, NULL, 'y'},
{"iterations", 1, NULL, 'n'},
{0, 0, NULL, 0}
@@ -329,10 +185,6 @@ int main(int argc, char *argv[])
geomfile = strdup(optarg);
break;
- case 'x' :
- prefix = strdup(optarg);
- break;
-
case 'j' :
nthreads = atoi(optarg);
break;
@@ -387,15 +239,6 @@ int main(int argc, char *argv[])
outfile = strdup("facetron.hkl");
}
- /* Sanitise prefix */
- if ( prefix == NULL ) {
- prefix = strdup("");
- } else {
- if ( config_checkprefix ) {
- prefix = check_prefix(prefix);
- }
- }
-
if ( sym == NULL ) sym = strdup("1");
/* Get detector geometry */
@@ -427,7 +270,11 @@ int main(int argc, char *argv[])
UnitCell *cell;
char *filename;
- char *fnamereal;
+ char *rval;
+ char line[1024];
+ RefList *peaks;
+ RefList *transfer;
+ Reflection *refl;
if ( find_chunk(fh, &cell, &filename) == 1 ) {
ERROR("Couldn't get all of the filenames and cells"
@@ -437,38 +284,83 @@ int main(int argc, char *argv[])
images[i].indexed_cell = cell;
- /* Mangle the filename now */
- if ( config_basename ) {
- char *tmp;
- tmp = safe_basename(filename);
- free(filename);
- filename = tmp;
- }
- fnamereal = malloc(1024);
- snprintf(fnamereal, 1023, "%s%s", prefix, filename);
- free(filename);
-
- images[i].filename = fnamereal;
+ images[i].filename = filename;
images[i].div = beam->divergence;
images[i].bw = beam->bandwidth;
images[i].det = det;
images[i].beam = beam;
images[i].osf = 1.0;
images[i].profile_radius = 0.005e9;
+ images[i].reflections = reflist_new();
/* Muppet proofing */
images[i].data = NULL;
images[i].flags = NULL;
- /* Get reflections from this image.
- * FIXME: Use the ones from the stream */
- integrate_image(&images[i], obs, sym);
+ /* Read integrated intensities from pattern */
+ peaks = reflist_new();
+ do {
+
+ Reflection *refl;
+ signed int h, k, l;
+ float intensity;
+ int r;
+
+ rval = fgets(line, 1023, fh);
+ chomp(line);
+
+ if ( (strlen(line) == 0) || (rval == NULL) ) break;
+
+ r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
+ if ( r != 4 ) continue;
+
+ /* Not interested in the central beam */
+ if ( (h==0) && (k==0) && (l==0) ) continue;
+
+ refl = add_refl(peaks, h, k, l);
+ set_int(refl, intensity);
+
+ } while ( (strlen(line) != 0) && (rval != NULL) );
+
+ /* Calculate initial partialities and fill in intensities from
+ * the stream */
+ transfer = find_intersections(&images[i], cell, 0);
+ images[i].reflections = reflist_new();
+
+ for ( refl = first_refl(transfer);
+ refl != NULL;
+ refl = next_refl(refl) ) {
+
+ Reflection *peak;
+ Reflection *new;
+ signed int h, k, l, ha, ka, la;
+ double r1, r2, p, x, y;
+ int clamp1, clamp2;
+
+ get_indices(refl, &h, &k, &l);
+ peak = find_refl(peaks, h, k, l);
+ if ( peak == NULL ) {
+ if ( (h==0) && (k==0) && (l==0) ) continue;
+ STATUS("%3i %3i %3i not found\n", h, k, l);
+ continue;
+ }
+
+ get_asymm(h, k, l, &ha, &ka, &la, sym);
+ add_item(obs, ha, ka, la);
+ new = add_refl(images[i].reflections, ha, ka, la);
+ get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
+ get_detector_pos(refl, &x, &y);
+ set_partial(new, r1, r2, p, clamp1, clamp2);
+ set_detector_pos(new, 0.0, x, y);
+
+ }
+ reflist_free(peaks);
+ reflist_free(transfer);
progress_bar(i, n_total_patterns-1, "Loading pattern data");
}
fclose(fh);
- free(prefix);
cts = new_list_count();
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 25664e69..93b65810 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -28,6 +28,28 @@
#include "cell.h"
+/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
+#define MAX_CYCLES (100)
+
+
+/* Refineable parameters */
+enum {
+ REF_SCALE,
+ REF_DIV,
+ REF_R,
+ REF_ASX,
+ REF_BSX,
+ REF_CSX,
+ REF_ASY,
+ REF_BSY,
+ REF_CSY,
+ REF_ASZ,
+ REF_BSZ,
+ REF_CSZ,
+ NUM_PARAMS
+};
+
+
/* Returns dp/dr at "r" */
static double partiality_gradient(double r, double profile_radius)
{
@@ -65,7 +87,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, Reflection *refl, double r)
+static double gradient(struct image *image, int k, Reflection *refl, double r)
{
double ds, tt, azi;
double nom, den;
@@ -187,7 +209,7 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift)
/* Apply the given shift to the 'k'th parameter of 'image'. */
-void apply_shift(struct image *image, int k, double shift)
+static void apply_shift(struct image *image, int k, double shift)
{
switch ( k ) {
@@ -226,67 +248,19 @@ void apply_shift(struct image *image, int k, double shift)
}
-double mean_partial_dev(struct image *image, RefList *reflections,
- const char *sym, double *i_full, FILE *graph)
-{
- int n_used;
- double delta_I = 0.0;
- Reflection *refl;
-
- n_used = 0;
- 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;
-
- get_indices(refl, &hind, &kind, &lind);
-
- if ( !get_scalable(refl) ) continue;
-
- /* 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;
- }
-
- get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
- I_full = lookup_intensity(i_full, ha, ka, la);
- 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/p, xc, yc,
- p, I_partial / I_full);
- }
-
- }
-
- return delta_I / (double)n_used;
-}
-
-
/* Perform one cycle of post refinement on 'image' against 'i_full' */
-double pr_iterate(struct image *image, double *i_full, const char *sym,
- RefList *reflections)
+static double pr_iterate(struct image *image, const double *i_full,
+ const char *sym)
{
gsl_matrix *M;
gsl_vector *v;
gsl_vector *shifts;
int param;
Reflection *refl;
+ RefList *reflections;
+ double max_shift;
+
+ reflections = image->reflections;
M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
v = gsl_vector_calloc(NUM_PARAMS);
@@ -355,14 +329,31 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
shifts = gsl_vector_alloc(NUM_PARAMS);
gsl_linalg_HH_solve(M, v, shifts);
+
+ max_shift = 0.0;
for ( param=0; param<NUM_PARAMS; param++ ) {
double shift = gsl_vector_get(shifts, param);
apply_shift(image, param, shift);
+ if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
}
gsl_matrix_free(M);
gsl_vector_free(v);
gsl_vector_free(shifts);
- return mean_partial_dev(image, reflections, sym, i_full, NULL);
+ return max_shift;
+}
+
+
+void pr_refine(struct image *image, const double *i_full, const char *sym)
+{
+ double max_shift;
+ int i;
+
+ i = 0;
+ do {
+ max_shift = pr_iterate(image, i_full, sym);
+ STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift);
+ i++;
+ } while ( (max_shift > 0.01) && (i < MAX_CYCLES) );
}
diff --git a/src/post-refinement.h b/src/post-refinement.h
index b3c6c074..a10fafe4 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -3,7 +3,7 @@
*
* Post refinement
*
- * (c) 2006-2010 Thomas White <taw@physics.org>
+ * (c) 2006-2011 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
@@ -21,36 +21,11 @@
#include <stdio.h>
#include "image.h"
-#include "reflist.h"
+#include "utils.h"
-/* Refineable parameters */
-enum {
- REF_SCALE,
- REF_DIV,
- REF_R,
- REF_ASX,
- REF_BSX,
- REF_CSX,
- REF_ASY,
- REF_BSY,
- REF_CSY,
- REF_ASZ,
- REF_BSZ,
- REF_CSZ,
- NUM_PARAMS
-};
-
-/* Apply the given shift to the 'k'th parameter of 'image'. */
-void apply_shift(struct image *image, int k, double shift);
-
-
-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,
- RefList *reflections);
+extern void pr_refine(struct image *image, const double *i_full,
+ const char *sym);
#endif /* POST_REFINEMENT_H */
diff --git a/src/reflist.c b/src/reflist.c
index d297c4d0..2bb4e5f6 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -19,7 +19,6 @@
struct _reflection {
/* Listy stuff */
- RefList *list;
unsigned int serial; /* Unique serial number, key */
struct _reflection *child[2]; /* Child nodes */
struct _reflection *parent; /* Parent node */
@@ -224,32 +223,6 @@ void set_partial(Reflection *refl, double r1, double r2, double p,
}
-void set_indices(Reflection *refl, signed int h, signed int k, signed int l)
-{
- /* Tricky, because the indices determine the position in the tree */
- Reflection copy;
- Reflection *new;
- Reflection transfer;
-
- /* Copy all data */
- copy = *refl;
-
- /* Delete and re-add with new indices */
- delete_refl(refl);
- new = add_refl(copy.list, h, k, l);
-
- /* Transfer data back */
- transfer = *new;
- *new = copy;
- new->list = transfer.list;
- new->parent = transfer.parent;
- new->child[0] = transfer.child[0];
- new->child[1] = transfer.child[1];
- new->h = transfer.h; new->k = transfer.k; new->l = transfer.l;
- new->serial = transfer.serial;
-}
-
-
void set_int(Reflection *refl, double intensity)
{
refl->intensity = intensity;
diff --git a/src/reflist.h b/src/reflist.h
index c7b2ad42..78aadf2a 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -46,8 +46,6 @@ 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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 13:44:48 +0000