aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/cell-utils.c19
-rw-r--r--libcrystfel/src/crystal.c64
-rw-r--r--libcrystfel/src/crystal.h12
-rw-r--r--libcrystfel/src/detector.c195
-rw-r--r--libcrystfel/src/detector.h8
-rw-r--r--libcrystfel/src/events.c12
-rw-r--r--libcrystfel/src/events.h12
-rw-r--r--libcrystfel/src/geometry.c339
-rw-r--r--libcrystfel/src/geometry.h34
-rw-r--r--libcrystfel/src/hdf5-file.c797
-rw-r--r--libcrystfel/src/hdf5-file.h21
-rw-r--r--libcrystfel/src/index.c34
-rw-r--r--libcrystfel/src/index.h8
-rw-r--r--libcrystfel/src/integration.c133
-rw-r--r--libcrystfel/src/integration.h7
-rw-r--r--libcrystfel/src/peaks.h4
-rw-r--r--libcrystfel/src/reflist-utils.c2
-rw-r--r--libcrystfel/src/reflist.c77
-rw-r--r--libcrystfel/src/reflist.h12
-rw-r--r--libcrystfel/src/render.c82
-rw-r--r--libcrystfel/src/render.h3
-rw-r--r--libcrystfel/src/stream.c57
-rw-r--r--libcrystfel/src/stream.h3
-rw-r--r--libcrystfel/src/utils.c157
-rw-r--r--libcrystfel/src/utils.h2
25 files changed, 1211 insertions, 883 deletions
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 5fc885e5..75419c99 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -3,13 +3,13 @@
*
* Unit Cell utility functions
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2009-2012,2014 Thomas White <taw@physics.org>
- * 2012 Lorenzo Galli
+ * 2009-2012,2014-2015 Thomas White <taw@physics.org>
+ * 2012 Lorenzo Galli
*
* This file is part of CrystFEL.
*
@@ -606,10 +606,12 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
/* "Un-center" the template unit cell to make the comparison easier */
template = uncenter_cell(template_in, &uncentering);
+ if ( template == NULL ) return NULL;
/* The candidate cell is also uncentered, because it might be centered
* if it came from (e.g.) MOSFLM */
cell = uncenter_cell(cell_in, NULL);
+ if ( cell == NULL ) return NULL;
if ( cell_get_reciprocal(template, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
@@ -1090,8 +1092,9 @@ static void determine_lattice(UnitCell *cell,
/**
* load_cell_from_pdb:
+ * @filename: The filename from which to load the cell
*
- * Loads a unit cell from a PDB file.
+ * Loads a unit cell from the CRYST1 line of a PDB file.
*
* Returns: a newly allocated %UnitCell.
*
@@ -1161,7 +1164,12 @@ UnitCell *load_cell_from_pdb(const char *filename)
fclose(fh);
- validate_cell(cell);
+ if ( cell != NULL ) {
+ validate_cell(cell);
+ } else {
+ ERROR("Failed to load cell from %s\n", filename);
+ }
+
return cell;
}
@@ -1225,6 +1233,7 @@ static int get_angle_rad(char **bits, int nbits, double *pl)
/**
* load_cell_from_file:
+ * @filename: The filename from which to load the cell
*
* Loads a unit cell from a file of any type (PDB or CrystFEL format)
*
diff --git a/libcrystfel/src/crystal.c b/libcrystfel/src/crystal.c
index c6e0a9dd..0c707331 100644
--- a/libcrystfel/src/crystal.c
+++ b/libcrystfel/src/crystal.c
@@ -3,11 +3,11 @@
*
* A class representing a single crystal
*
- * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2013-2015 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2013 Thomas White <taw@physics.org>
+ * 2013-2015 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -56,6 +56,7 @@ struct _crystal
UnitCell *cell;
double m; /* Mosaicity in radians */
double osf;
+ double Bfac;
double profile_radius;
int pr_dud;
double resolution_limit;
@@ -68,6 +69,9 @@ struct _crystal
/* User flag, e.g. for "this is a bad crystal". */
int user_flag;
+
+ /* Text notes, which go in the stream */
+ char *notes;
};
@@ -91,9 +95,11 @@ Crystal *crystal_new()
cryst->cell = NULL;
cryst->reflections = NULL;
- cryst->resolution_limit = 0.0;
+ cryst->resolution_limit = INFINITY;
cryst->n_saturated = 0;
cryst->n_implausible = 0;
+ cryst->notes = NULL;
+ cryst->user_flag = 0;
return cryst;
}
@@ -188,6 +194,12 @@ double crystal_get_osf(Crystal *cryst)
}
+double crystal_get_Bfac(Crystal *cryst)
+{
+ return cryst->Bfac;
+}
+
+
int crystal_get_user_flag(Crystal *cryst)
{
return cryst->user_flag;
@@ -200,6 +212,12 @@ double crystal_get_mosaicity(Crystal *cryst)
}
+const char *crystal_get_notes(Crystal *cryst)
+{
+ return cryst->notes;
+}
+
+
/********************************** Setters ***********************************/
@@ -251,6 +269,12 @@ void crystal_set_osf(Crystal *cryst, double osf)
}
+void crystal_set_Bfac(Crystal *cryst, double Bfac)
+{
+ cryst->Bfac = Bfac;
+}
+
+
void crystal_set_user_flag(Crystal *cryst, int user_flag)
{
cryst->user_flag = user_flag;
@@ -261,3 +285,35 @@ void crystal_set_mosaicity(Crystal *cryst, double m)
{
cryst->m = m;
}
+
+
+void crystal_set_notes(Crystal *cryst, const char *notes)
+{
+ free(cryst->notes); /* free(NULL) is OK */
+ cryst->notes = strdup(notes);
+}
+
+
+void crystal_add_notes(Crystal *cryst, const char *notes_add)
+{
+ size_t len;
+ char *nnotes;
+
+ if ( cryst->notes == NULL ) {
+ crystal_set_notes(cryst, notes_add);
+ return;
+ }
+
+ len = strlen(notes_add) + strlen(cryst->notes) + 2;
+ nnotes = malloc(len);
+ if ( nnotes == NULL ) {
+ ERROR("Failed to add notes to crystal.\n");
+ return;
+ }
+
+ strcpy(nnotes, cryst->notes);
+ strcat(nnotes, "\n");
+ strcat(nnotes, notes_add);
+ free(cryst->notes);
+ cryst->notes = nnotes;
+}
diff --git a/libcrystfel/src/crystal.h b/libcrystfel/src/crystal.h
index 987cabbf..f0d732ff 100644
--- a/libcrystfel/src/crystal.h
+++ b/libcrystfel/src/crystal.h
@@ -3,11 +3,11 @@
*
* A class representing a single crystal
*
- * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2013-2015 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2013 Thomas White <taw@physics.org>
+ * 2013-2015 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -62,8 +62,10 @@ extern long long int crystal_get_num_saturated_reflections(Crystal *cryst);
extern long long int crystal_get_num_implausible_reflections(Crystal *cryst);
extern int crystal_get_user_flag(Crystal *cryst);
extern double crystal_get_osf(Crystal *cryst);
+extern double crystal_get_Bfac(Crystal *cryst);
extern struct image *crystal_get_image(Crystal *cryst);
extern double crystal_get_mosaicity(Crystal *cryst);
+extern const char *crystal_get_notes(Crystal *cryst);
extern void crystal_set_cell(Crystal *cryst, UnitCell *cell);
extern void crystal_set_profile_radius(Crystal *cryst, double r);
@@ -75,8 +77,12 @@ extern void crystal_set_num_implausible_reflections(Crystal *cryst,
long long int n);
extern void crystal_set_user_flag(Crystal *cryst, int flag);
extern void crystal_set_osf(Crystal *cryst, double osf);
+extern void crystal_set_Bfac(Crystal *cryst, double B);
extern void crystal_set_image(Crystal *cryst, struct image *image);
extern void crystal_set_mosaicity(Crystal *cryst, double m);
+extern void crystal_set_notes(Crystal *cryst, const char *notes);
+
+extern void crystal_add_notes(Crystal *cryst, const char *notes_add);
#ifdef __cplusplus
}
diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
index a0de3be9..7979315b 100644
--- a/libcrystfel/src/detector.c
+++ b/libcrystfel/src/detector.c
@@ -3,12 +3,12 @@
*
* Detector properties
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
* Authors:
- * 2009-2014 Thomas White <taw@physics.org>
+ * 2009-2015 Thomas White <taw@physics.org>
* 2014 Valerio Mariani
* 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
* 2011 Andrew Aquila
@@ -339,16 +339,30 @@ int in_bad_region(struct detector *det, double fs, double ss)
struct badregion *b = &det->bad[i];
+ if ( (b->panel != NULL)
+ && (strcmp(b->panel, p->name) != 0) ) continue;
+
if ( b->is_fsss ) {
- if ( fs < b->min_fs ) continue;
- if ( fs > b->max_fs ) continue;
- if ( ss < b->min_ss ) continue;
- if ( ss > b->max_ss ) continue;
+
+ int nfs, nss;
+
+ /* fs/ss bad regions are specified according to the
+ * original coordinates */
+ nfs = (fs-p->min_fs) + p->orig_min_fs;
+ nss = (ss-p->min_ss) + p->orig_min_ss;
+
+ if ( nfs < b->min_fs ) continue;
+ if ( nfs > b->max_fs ) continue;
+ if ( nss < b->min_ss ) continue;
+ if ( nss > b->max_ss ) continue;
+
} else {
+
if ( rx < b->min_x ) continue;
if ( rx > b->max_x ) continue;
if ( ry < b->min_y ) continue;
if ( ry > b->max_y ) continue;
+
}
return 1;
@@ -386,6 +400,18 @@ double get_tt(struct image *image, double fs, double ss, int *err)
}
+int detector_has_clen_references(struct detector *det)
+{
+ int i;
+
+ for ( i=0; i<det->n_panels; i++ ) {
+ if ( det->panels[i].clen_from != NULL ) return 1;
+ }
+
+ return 0;
+}
+
+
void record_image(struct image *image, int do_poisson, double background,
gsl_rng *rng, double beam_radius, double nphotons)
{
@@ -558,13 +584,17 @@ void fill_in_values(struct detector *det, struct hdfile *f, struct event* ev)
if ( p->clen_from != NULL ) {
+ double val;
+ int r;
- if (det->path_dim !=0 || det->dim_dim !=0 ){
- p->clen = get_ev_based_value(f, p->clen_from,
- ev) * 1.0e-3;
+ r = hdfile_get_value(f, p->clen_from, ev, &val,
+ H5T_NATIVE_DOUBLE);
+ if ( r ) {
+ ERROR("Failed to read '%s'\n", p->clen_from);
} else {
- p->clen = get_value(f, p->clen_from) * 1.0e-3;
+ p->clen = val * 1.0e-3;
}
+
}
p->clen += p->coffset;
@@ -681,7 +711,8 @@ static struct badregion *new_bad_region(struct detector *det, const char *name)
new->max_fs = 0;
new->min_ss = 0;
new->max_ss = 0;
- new->is_fsss = 0;
+ new->is_fsss = 99; /* Slightly nasty: means "unassigned" */
+ new->panel = NULL;
strcpy(new->name, name);
return new;
@@ -952,43 +983,54 @@ static int parse_field_for_panel(struct panel *panel, const char *key,
}
-static int parse_field_bad(struct badregion *panel, const char *key,
- const char *val)
+static int check_badr_fsss(struct badregion *badr, int is_fsss)
+{
+ /* First assignment? */
+ if ( badr->is_fsss == 99 ) {
+ badr->is_fsss = is_fsss;
+ return 0;
+ }
+
+ if ( is_fsss != badr->is_fsss ) {
+ ERROR("You can't mix x/y and fs/ss in a bad region.\n");
+ return 1;
+ }
+
+ return 0;
+}
+
+
+static int parse_field_bad(struct badregion *badr, const char *key,
+ const char *val)
{
int reject = 0;
if ( strcmp(key, "min_x") == 0 ) {
- panel->min_x = atof(val);
- if ( panel->is_fsss ) {
- ERROR("You can't mix x/y and fs/ss in a bad region.\n");
- }
+ badr->min_x = atof(val);
+ reject = check_badr_fsss(badr, 0);
} else if ( strcmp(key, "max_x") == 0 ) {
- panel->max_x = atof(val);
- if ( panel->is_fsss ) {
- ERROR("You can't mix x/y and fs/ss in a bad region.\n");
- }
+ badr->max_x = atof(val);
+ reject = check_badr_fsss(badr, 0);
} else if ( strcmp(key, "min_y") == 0 ) {
- panel->min_y = atof(val);
- if ( panel->is_fsss ) {
- ERROR("You can't mix x/y and fs/ss in a bad region.\n");
- }
+ badr->min_y = atof(val);
+ reject = check_badr_fsss(badr, 0);
} else if ( strcmp(key, "max_y") == 0 ) {
- panel->max_y = atof(val);
- if ( panel->is_fsss ) {
- ERROR("You can't mix x/y and fs/ss in a bad region.\n");
- }
+ badr->max_y = atof(val);
+ reject = check_badr_fsss(badr, 0);
} else if ( strcmp(key, "min_fs") == 0 ) {
- panel->min_fs = atof(val);
- panel->is_fsss = 1;
+ badr->min_fs = atof(val);
+ reject = check_badr_fsss(badr, 1);
} else if ( strcmp(key, "max_fs") == 0 ) {
- panel->max_fs = atof(val);
- panel->is_fsss = 1;
+ badr->max_fs = atof(val);
+ reject = check_badr_fsss(badr, 1);
} else if ( strcmp(key, "min_ss") == 0 ) {
- panel->min_ss = atof(val);
- panel->is_fsss = 1;
+ badr->min_ss = atof(val);
+ reject = check_badr_fsss(badr, 1);
} else if ( strcmp(key, "max_ss") == 0 ) {
- panel->max_ss = atof(val);
- panel->is_fsss = 1;
+ badr->max_ss = atof(val);
+ reject = check_badr_fsss(badr, 1);
+ } else if ( strcmp(key, "panel") == 0 ) {
+ badr->panel = strdup(val);
} else {
ERROR("Unrecognised field '%s'\n", key);
}
@@ -1026,12 +1068,14 @@ static void parse_toplevel(struct detector *det, struct beam_params *beam,
det->defaults.coffset = atof(val);
} else if ( strcmp(key, "photon_energy") == 0 ) {
- if ( strncmp(val, "/", 1) == 0 ) {
- beam->photon_energy = 0.0;
- beam->photon_energy_from = strdup(val);
- } else {
- beam->photon_energy = atof(val);
- beam->photon_energy_from = NULL;
+ if ( beam != NULL ) {
+ if ( strncmp(val, "/", 1) == 0 ) {
+ beam->photon_energy = 0.0;
+ beam->photon_energy_from = strdup(val);
+ } else {
+ beam->photon_energy = atof(val);
+ beam->photon_energy_from = NULL;
+ }
}
} else if ( strcmp(key, "photon_energy_scale") == 0 ) {
@@ -1191,7 +1235,7 @@ struct detector *get_detector_geometry(const char *filename,
det->defaults.orig_max_ss = -1;
det->defaults.cnx = NAN;
det->defaults.cny = NAN;
- det->defaults.clen = -1.0;
+ det->defaults.clen = NAN;
det->defaults.coffset = 0.0;
det->defaults.res = -1.0;
det->defaults.badrow = '-';
@@ -1464,7 +1508,7 @@ struct detector *get_detector_geometry(const char *filename,
" panel %s\n", det->panels[i].name);
reject = 1;
}
- if ( (det->panels[i].clen < 0.0)
+ if ( isnan(det->panels[i].clen)
&& (det->panels[i].clen_from == NULL) ) {
ERROR("Please specify the camera length for"
" panel %s\n", det->panels[i].name);
@@ -1511,24 +1555,8 @@ struct detector *get_detector_geometry(const char *filename,
}
for ( i=0; i<det->n_bad; i++ ) {
-
- if ( !det->bad[i].is_fsss && isnan(det->bad[i].min_x) ) {
- ERROR("Please specify the minimum x coordinate for"
- " bad region %s\n", det->bad[i].name);
- reject = 1;
- }
- if ( !det->bad[i].is_fsss && isnan(det->bad[i].min_y) ) {
- ERROR("Please specify the minimum y coordinate for"
- " bad region %s\n", det->bad[i].name);
- reject = 1;
- }
- if ( !det->bad[i].is_fsss && isnan(det->bad[i].max_x) ) {
- ERROR("Please specify the maximum x coordinate for"
- " bad region %s\n", det->bad[i].name);
- reject = 1;
- }
- if ( !det->bad[i].is_fsss && isnan(det->bad[i].max_y) ) {
- ERROR("Please specify the maximum y coordinate for"
+ if ( det->bad[i].is_fsss == 99 ) {
+ ERROR("Please specify the coordinate ranges for"
" bad region %s\n", det->bad[i].name);
reject = 1;
}
@@ -1569,6 +1597,7 @@ out:
if ( p == NULL ) {
ERROR("Cannot add panel to rigid group\n");
ERROR("Panel not found: %s\n", bits[pi]);
+ return NULL;
}
add_to_rigid_group(rigidgroup, p);
free(bits[pi]);
@@ -1596,7 +1625,8 @@ out:
r = find_rigid_group_by_name(det, bits[rgi]);
if ( r == NULL ) {
ERROR("Cannot add rigid group to collection\n");
- ERROR("Rigid groups not found: %s\n", bits[rgi]);
+ ERROR("Rigid group not found: %s\n", bits[rgi]);
+ return NULL;
}
add_to_rigid_group_coll(rgcollection, r);
free(bits[rgi]);
@@ -1703,6 +1733,8 @@ struct detector *copy_geom(const struct detector *in)
out->n_rigid_groups = 0;
out->rigid_groups = NULL;
+ out->n_rg_collections = 0;
+ out->rigid_group_collections = NULL;
for ( i=0; i<out->n_panels; i++ ) {
@@ -1719,7 +1751,7 @@ struct detector *copy_geom(const struct detector *in)
if ( p->data != NULL ) {
/* Make a copy of the data fields unique to this
* copy of the structure. */
- p->clen_from = strdup(p->clen_from);
+ p->data = strdup(p->data);
}
if ( p->dim_structure != NULL ) {
@@ -1753,14 +1785,25 @@ struct detector *copy_geom(const struct detector *in)
int rgi;
for ( rgi=0; rgi<in->n_rigid_groups; rgi++ ) {
+
if ( panel_is_in_rigid_group(in->rigid_groups[rgi],
- &in->panels[i])) {
- add_to_rigid_group(find_or_add_rg(out,
- in->rigid_groups[rgi]->name),
- &out->panels[i]);
+ &in->panels[i]) )
+ {
+ struct rigid_group *g;
+ g = find_or_add_rg(out,
+ in->rigid_groups[rgi]->name);
+ add_to_rigid_group(g, &out->panels[i]);
}
}
+ if ( &in->panels[i] == in->furthest_out_panel ) {
+ out->furthest_out_panel = &out->panels[i];
+ }
+ if ( &in->panels[i] == in->furthest_in_panel ) {
+ out->furthest_in_panel = &out->panels[i];
+ }
+
+
}
for ( i=0; i<in->n_rigid_groups; i++ ) {
@@ -1768,13 +1811,17 @@ struct detector *copy_geom(const struct detector *in)
int rgci;
for ( rgci=0; rgci<in->n_rg_collections; rgci++ ) {
+
+ const char *n = in->rigid_group_collections[rgci]->name;
+
if ( rigid_group_is_in_collection(
in->rigid_group_collections[rgci],
in->rigid_groups[i]) )
{
- add_to_rigid_group_coll(find_or_add_rg_coll(out,
- in->rigid_group_collections[rgci]->name),
- out->rigid_groups[i]);
+ struct rg_collection *rgcoll;
+ rgcoll = find_or_add_rg_coll(out, n);
+ add_to_rigid_group_coll(rgcoll,
+ out->rigid_groups[i]);
}
}
}
@@ -1894,7 +1941,8 @@ static void process_panel_fields(const struct panel *p, char *line,
if(strstr(bits[1], "fs") != NULL &&
strstr(bits[1], "min_fs") == NULL &&
- strstr(bits[1], "max_fs") == NULL) {
+ strstr(bits[1], "max_fs") == NULL &&
+ strstr(bits[1], "offset") == NULL ) {
sprintf(string_to_write, "%+fx %+fy",
p->fsx, p->fsy);
@@ -2030,7 +2078,6 @@ int write_detector_geometry_2(const char *geometry_filename,
fh = fopen(output_filename, "w");
if ( fh == NULL ) {
- fclose(fh);
return 1;
}
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
index a82134bf..acb6609f 100644
--- a/libcrystfel/src/detector.h
+++ b/libcrystfel/src/detector.h
@@ -3,12 +3,12 @@
*
* Detector properties
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
* Authors:
- * 2009-2014 Thomas White <taw@physics.org>
+ * 2009-2015 Thomas White <taw@physics.org>
* 2011-2012 Richard Kirian <rkirian@asu.edu>
* 2014 Valerio Mariani
* 2011 Andrew Aquila
@@ -38,6 +38,7 @@
#define DETECTOR_H
struct rigid_group;
+struct rg_collection;
struct detector;
struct panel;
struct badregion;
@@ -131,6 +132,7 @@ struct badregion
{
char name[1024];
int is_fsss;
+ char *panel;
double min_x;
double max_x;
@@ -253,6 +255,8 @@ extern int single_panel_data_source (struct detector *det, const char *element);
struct rg_collection *find_rigid_group_collection_by_name(struct detector *det,
const char *name);
+extern int detector_has_clen_references(struct detector *det);
+
#ifdef __cplusplus
}
#endif
diff --git a/libcrystfel/src/events.c b/libcrystfel/src/events.c
index b008eee6..731f39f8 100644
--- a/libcrystfel/src/events.c
+++ b/libcrystfel/src/events.c
@@ -292,7 +292,7 @@ char *get_event_string(struct event *ev)
char *new_ret_string;
int ret_string_len;
- if ( ev == NULL ) return "(none)";
+ if ( ev == NULL ) return strdup("(none)");
if ( ev->path_length != 0 ) {
@@ -356,7 +356,7 @@ char *get_event_string(struct event *ev)
strncpy(&ret_string[ret_string_len],"/", 1);
strncpy(&ret_string[ret_string_len+1], num_buf,
- strlen(num_buf));
+ strlen(num_buf));
ret_string_len += 1+strlen(num_buf);
}
@@ -412,16 +412,13 @@ struct event *get_event_from_event_string(const char *ev_string)
if ( ev == NULL ) return NULL;
if ( strlen(buf_path) !=0 ) {
-
+ start = buf_path;
do {
-
- start = buf_path;
-
char buf[2014];
sep = strstr(start, "/");
- if ( sep != NULL ) {
+ if ( sep != NULL ) {
strncpy(buf, start, sep-start);
buf[sep-start]='\0';
push_path_entry_to_event(ev, buf);
@@ -462,6 +459,7 @@ struct event *get_event_from_event_string(const char *ev_string)
push_dim_entry_to_event(ev, buf_int);
}
+
} while (sep);
}
diff --git a/libcrystfel/src/events.h b/libcrystfel/src/events.h
index 8cb00962..7f9c6731 100644
--- a/libcrystfel/src/events.h
+++ b/libcrystfel/src/events.h
@@ -67,7 +67,7 @@ struct dim_structure
int num_dims;
};
-extern struct event *initialize_event();
+extern struct event *initialize_event(void);
extern int push_path_entry_to_event(struct event *ev, const char *entry);
extern int pop_path_entry_from_event(struct event *ev);
extern int push_dim_entry_to_event(struct event *ev, int entry);
@@ -82,22 +82,22 @@ extern char *partial_event_substitution(struct event *ev, const char *data);
extern char *retrieve_full_path(struct event *ev, const char *data);
-extern struct filename_plus_event *initialize_filename_plus_event();
+extern struct filename_plus_event *initialize_filename_plus_event(void);
extern void free_filename_plus_event(struct filename_plus_event *fpe);
-extern struct event_list *initialize_event_list();
+extern struct event_list *initialize_event_list(void);
extern int append_event_to_event_list(struct event_list *ev_list,
struct event *ev);
-int add_non_existing_event_to_event_list(struct event_list *ev_list,
+extern int add_non_existing_event_to_event_list(struct event_list *ev_list,
struct event *ev);
extern struct event_list *copy_event_list(struct event_list *el);
extern int find_event(struct event *ev, struct event_list *el);
extern void free_event_list(struct event_list *el);
-extern struct dim_structure *initialize_dim_structure();
-extern struct dim_structure *default_dim_structure();
+extern struct dim_structure *initialize_dim_structure(void);
+extern struct dim_structure *default_dim_structure(void);
extern int set_dim_structure_entry(struct dim_structure *hsd,
const char *string_dim,
const char *val_string);
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index b55b696d..acfc08ee 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -47,58 +47,65 @@
#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;
- p = &det->panels[i];
+ 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;
- /* Coordinates of peak relative to central beam, in m */
- xd = p->clen * x / den;
- yd = p->clen * y / den;
+ /* 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;
- /* Convert to pixels */
- xd *= p->res;
- yd *= p->res;
+ return 1;
+}
- /* Convert to relative to the panel corner */
- plx = xd - p->cnx;
- ply = yd - p->cny;
+static signed int locate_peak(double x, double y, double z, double k,
+ struct detector *det, double *pfs, double *pss)
+{
+ int i;
+
+ *pfs = -1; *pss = -1;
+
+ for ( i=0; i<det->n_panels; i++ ) {
- fs = p->xfs*plx + p->yfs*ply;
- ss = p->xss*plx + p->yss*ply;
+ struct panel *p;
- fs += p->min_fs;
- ss += p->min_ss;
+ 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 +189,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 +203,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 +215,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 +228,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 ( (image->det != NULL) && (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);
+
+ }
+
+ /* 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) ) {
- /* 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);
+ 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 +276,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);
@@ -266,6 +296,92 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
}
+double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
+{
+ double azi;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double xl, yl, zl;
+ signed int hs, ks, ls;
+ double rlow, rhigh, p;
+ double philow, phihigh, phi;
+ double khigh, klow;
+ double tl, cet, cez;
+
+ get_partial(refl, &rlow, &rhigh, &p);
+
+ get_symmetric_indices(refl, &hs, &ks, &ls);
+
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ xl = hs*asx + ks*bsx + ls*csx;
+ yl = hs*asy + ks*bsy + ls*csy;
+ zl = hs*asz + ks*bsz + ls*csz;
+
+ /* "low" gives the largest Ewald sphere (wavelength short => k large)
+ * "high" gives the smallest Ewald sphere (wavelength long => k small)
+ */
+ klow = 1.0/(image->lambda - image->lambda*image->bw/2.0);
+ khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0);
+
+ tl = sqrt(xl*xl + yl*yl);
+
+ cet = -sin(image->div/2.0) * klow;
+ cez = -cos(image->div/2.0) * klow;
+ philow = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
+
+ cet = -sin(image->div/2.0) * khigh;
+ cez = -cos(image->div/2.0) * khigh;
+ phihigh = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
+
+ /* Approximation: philow and phihigh are very similar */
+ phi = (philow + phihigh) / 2.0;
+
+ azi = atan2(yl, xl);
+
+ switch ( k ) {
+
+ case GPARAM_ASX :
+ return - hs * sin(phi) * cos(azi);
+
+ case GPARAM_BSX :
+ return - ks * sin(phi) * cos(azi);
+
+ case GPARAM_CSX :
+ return - ls * sin(phi) * cos(azi);
+
+ case GPARAM_ASY :
+ return - hs * sin(phi) * sin(azi);
+
+ case GPARAM_BSY :
+ return - ks * sin(phi) * sin(azi);
+
+ case GPARAM_CSY :
+ return - ls * sin(phi) * sin(azi);
+
+ case GPARAM_ASZ :
+ return - hs * cos(phi);
+
+ case GPARAM_BSZ :
+ return - ks * cos(phi);
+
+ case GPARAM_CSZ :
+ return - ls * cos(phi);
+
+ case GPARAM_DETX :
+ case GPARAM_DETY :
+ case GPARAM_CLEN :
+ return 0.0;
+
+ }
+
+ ERROR("No r gradient defined for parameter %i\n", k);
+ abort();
+}
+
+
RefList *find_intersections(struct image *image, Crystal *cryst,
PartialityModel pmodel)
{
@@ -339,7 +455,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 +469,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 +485,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 +493,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 +506,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..152c0e47 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.
@@ -61,19 +61,37 @@ typedef enum {
} PartialityModel;
+
+/* Enumeration of parameters which may want to be refined */
+enum gparam {
+ GPARAM_ASX,
+ GPARAM_ASY,
+ GPARAM_ASZ,
+ GPARAM_BSX,
+ GPARAM_BSY,
+ GPARAM_BSZ,
+ GPARAM_CSX,
+ GPARAM_CSY,
+ GPARAM_CSZ,
+ GPARAM_R,
+ GPARAM_DIV,
+ GPARAM_DETX,
+ GPARAM_DETY,
+ GPARAM_CLEN,
+ GPARAM_OSF,
+ GPARAM_BFAC
+};
+
+
extern RefList *find_intersections(struct image *image, Crystal *cryst,
PartialityModel pmodel);
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 double r_gradient(UnitCell *cell, int k, Reflection *refl,
+ struct image *image);
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/hdf5-file.c b/libcrystfel/src/hdf5-file.c
index 85d4af72..3a08884a 100644
--- a/libcrystfel/src/hdf5-file.c
+++ b/libcrystfel/src/hdf5-file.c
@@ -378,169 +378,170 @@ static float *read_hdf5_data(struct hdfile *f, char *path, int line)
}
-int get_peaks(struct image *image, struct hdfile *f, const char *p,
- int cxi_format, struct filename_plus_event *fpe)
+/* Get peaks from HDF5, in "CXI format" (as in "CXIDB") */
+int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
+ struct filename_plus_event *fpe)
{
- if ( cxi_format ) {
-
- char path_n[1024];
- char path_x[1024];
- char path_y[1024];
- char path_i[1024];
- int r;
- int pk;
+ char path_n[1024];
+ char path_x[1024];
+ char path_y[1024];
+ char path_i[1024];
+ int r;
+ int pk;
- int line = 0;
- int num_peaks;
+ int line = 0;
+ int num_peaks;
- float *buf_x;
- float *buf_y;
- float *buf_i;
+ float *buf_x;
+ float *buf_y;
+ float *buf_i;
- if ( (fpe != NULL) && (fpe->ev != NULL)
- && (fpe->ev->dim_entries != NULL) )
- {
- line = fpe->ev->dim_entries[0];
- } else {
- ERROR("CXI format peak list format selected,"
- "but file has no event structure");
- return 1;
- }
-
- snprintf(path_n, 1024, "%s/nPeaks", p);
- snprintf(path_x, 1024, "%s/peakXPosRaw", p);
- snprintf(path_y, 1024, "%s/peakYPosRaw", p);
- snprintf(path_i, 1024, "%s/peakTotalIntensity", p);
+ if ( (fpe != NULL) && (fpe->ev != NULL)
+ && (fpe->ev->dim_entries != NULL) )
+ {
+ line = fpe->ev->dim_entries[0];
+ } else {
+ ERROR("CXI format peak list format selected,"
+ "but file has no event structure");
+ return 1;
+ }
- r = read_peak_count(f, path_n, line, &num_peaks);
- if ( r != 0 ) return 1;
+ snprintf(path_n, 1024, "%s/nPeaks", p);
+ snprintf(path_x, 1024, "%s/peakXPosRaw", p);
+ snprintf(path_y, 1024, "%s/peakYPosRaw", p);
+ snprintf(path_i, 1024, "%s/peakTotalIntensity", p);
- buf_x = read_hdf5_data(f, path_x, line);
- if ( r != 0 ) return 1;
+ r = read_peak_count(f, path_n, line, &num_peaks);
+ if ( r != 0 ) return 1;
- buf_y = read_hdf5_data(f, path_y, line);
- if ( r != 0 ) return 1;
+ buf_x = read_hdf5_data(f, path_x, line);
+ if ( r != 0 ) return 1;
- buf_i = read_hdf5_data(f, path_i, line);
- if ( r != 0 ) return 1;
+ buf_y = read_hdf5_data(f, path_y, line);
+ if ( r != 0 ) return 1;
- if ( image->features != NULL ) {
- image_feature_list_free(image->features);
- }
- image->features = image_feature_list_new();
+ buf_i = read_hdf5_data(f, path_i, line);
+ if ( r != 0 ) return 1;
- for ( pk=0; pk<num_peaks; pk++ ) {
+ if ( image->features != NULL ) {
+ image_feature_list_free(image->features);
+ }
+ image->features = image_feature_list_new();
- float fs, ss, val;
- struct panel *p;
+ for ( pk=0; pk<num_peaks; pk++ ) {
- fs = buf_x[pk];
- ss = buf_y[pk];
- val = buf_i[pk];
+ float fs, ss, val;
+ struct panel *p;
- p = find_orig_panel(image->det, fs, ss);
- if ( p == NULL ) continue;
- if ( p->no_index ) continue;
+ fs = buf_x[pk];
+ ss = buf_y[pk];
+ val = buf_i[pk];
- /* Convert coordinates to match rearranged
- * panels in memory */
- fs = fs - p->orig_min_fs + p->min_fs;
- ss = ss - p->orig_min_ss + p->min_ss;
+ p = find_orig_panel(image->det, fs, ss);
+ if ( p == NULL ) continue;
+ if ( p->no_index ) continue;
- image_add_feature(image->features, fs, ss, image,
- val, NULL);
+ /* Convert coordinates to match rearranged
+ * panels in memory */
+ fs = fs - p->orig_min_fs + p->min_fs;
+ ss = ss - p->orig_min_ss + p->min_ss;
- }
+ image_add_feature(image->features, fs, ss, image,
+ val, NULL);
- } else {
+ }
- hid_t dh, sh;
- hsize_t size[2];
- hsize_t max_size[2];
- int i;
- float *buf;
- herr_t r;
- int tw;
+ return 0;
+}
- dh = H5Dopen2(f->fh, p, H5P_DEFAULT);
- if ( dh < 0 ) {
- ERROR("Peak list (%s) not found.\n", p);
- return 1;
- }
- sh = H5Dget_space(dh);
- if ( sh < 0 ) {
- H5Dclose(dh);
- ERROR("Couldn't get dataspace for peak list.\n");
- return 1;
- }
+int get_peaks(struct image *image, struct hdfile *f, const char *p)
+{
+ hid_t dh, sh;
+ hsize_t size[2];
+ hsize_t max_size[2];
+ int i;
+ float *buf;
+ herr_t r;
+ int tw;
- if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
- ERROR("Peak list has the wrong dimensionality (%i).\n",
- H5Sget_simple_extent_ndims(sh));
- H5Sclose(sh);
- H5Dclose(dh);
- return 1;
- }
+ dh = H5Dopen2(f->fh, p, H5P_DEFAULT);
+ if ( dh < 0 ) {
+ ERROR("Peak list (%s) not found.\n", p);
+ return 1;
+ }
- H5Sget_simple_extent_dims(sh, size, max_size);
+ sh = H5Dget_space(dh);
+ if ( sh < 0 ) {
+ H5Dclose(dh);
+ ERROR("Couldn't get dataspace for peak list.\n");
+ return 1;
+ }
- tw = size[1];
- if ( (tw != 3) && (tw != 4) ) {
- H5Sclose(sh);
- H5Dclose(dh);
- ERROR("Peak list has the wrong dimensions.\n");
- return 1;
- }
+ if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
+ ERROR("Peak list has the wrong dimensionality (%i).\n",
+ H5Sget_simple_extent_ndims(sh));
+ H5Sclose(sh);
+ H5Dclose(dh);
+ return 1;
+ }
- buf = malloc(sizeof(float)*size[0]*size[1]);
- if ( buf == NULL ) {
- H5Sclose(sh);
- H5Dclose(dh);
- ERROR("Couldn't reserve memory for the peak list.\n");
- return 1;
- }
- r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, buf);
- if ( r < 0 ) {
- ERROR("Couldn't read peak list.\n");
- free(buf);
- return 1;
- }
+ H5Sget_simple_extent_dims(sh, size, max_size);
- if ( image->features != NULL ) {
- image_feature_list_free(image->features);
- }
- image->features = image_feature_list_new();
+ tw = size[1];
+ if ( (tw != 3) && (tw != 4) ) {
+ H5Sclose(sh);
+ H5Dclose(dh);
+ ERROR("Peak list has the wrong dimensions.\n");
+ return 1;
+ }
- for ( i=0; i<size[0]; i++ ) {
+ buf = malloc(sizeof(float)*size[0]*size[1]);
+ if ( buf == NULL ) {
+ H5Sclose(sh);
+ H5Dclose(dh);
+ ERROR("Couldn't reserve memory for the peak list.\n");
+ return 1;
+ }
+ r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, buf);
+ if ( r < 0 ) {
+ ERROR("Couldn't read peak list.\n");
+ free(buf);
+ return 1;
+ }
- float fs, ss, val;
- struct panel *p;
+ if ( image->features != NULL ) {
+ image_feature_list_free(image->features);
+ }
+ image->features = image_feature_list_new();
- fs = buf[tw*i+0];
- ss = buf[tw*i+1];
- val = buf[tw*i+2];
+ for ( i=0; i<size[0]; i++ ) {
- p = find_orig_panel(image->det, fs, ss);
- if ( p == NULL ) continue;
- if ( p->no_index ) continue;
+ float fs, ss, val;
+ struct panel *p;
- /* Convert coordinates to match rearranged panels in memory */
- fs = fs - p->orig_min_fs + p->min_fs;
- ss = ss - p->orig_min_ss + p->min_ss;
+ fs = buf[tw*i+0];
+ ss = buf[tw*i+1];
+ val = buf[tw*i+2];
- image_add_feature(image->features, fs, ss, image, val,
- NULL);
+ p = find_orig_panel(image->det, fs, ss);
+ if ( p == NULL ) continue;
+ if ( p->no_index ) continue;
- }
+ /* Convert coordinates to match rearranged panels in memory */
+ fs = fs - p->orig_min_fs + p->min_fs;
+ ss = ss - p->orig_min_ss + p->min_ss;
- free(buf);
- H5Sclose(sh);
- H5Dclose(dh);
+ image_add_feature(image->features, fs, ss, image, val,
+ NULL);
}
+ free(buf);
+ H5Sclose(sh);
+ H5Dclose(dh);
+
return 0;
}
@@ -1150,6 +1151,201 @@ static int unpack_panels(struct image *image, struct detector *det)
}
+static int get_scalar_value(struct hdfile *f, const char *name, void *val,
+ hid_t memtype)
+{
+ hid_t dh;
+ hid_t type;
+ hid_t class;
+ herr_t r;
+ int check;
+
+ if ( !hdfile_is_scalar(f, name, 1) ) return 1;
+
+ check = check_path_existence(f->fh, name);
+ if ( check == 0 ) {
+ ERROR("No such float field '%s'\n", name);
+ return 1;
+ }
+
+ dh = H5Dopen2(f->fh, name, H5P_DEFAULT);
+
+ type = H5Dget_type(dh);
+ class = H5Tget_class(type);
+
+ if ( class != H5T_FLOAT ) {
+ ERROR("Not a floating point value.\n");
+ H5Tclose(type);
+ H5Dclose(dh);
+ return 1;
+ }
+
+ r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, val);
+ if ( r < 0 ) {
+ ERROR("Couldn't read value.\n");
+ H5Tclose(type);
+ H5Dclose(dh);
+ return 1;
+ }
+
+ return 0;
+}
+
+
+static int get_ev_based_value(struct hdfile *f, const char *name,
+ struct event *ev, void *val, hid_t memtype)
+{
+ hid_t dh;
+ hid_t type;
+ hid_t class;
+ hid_t sh;
+ hid_t ms;
+ hsize_t *f_offset = NULL;
+ hsize_t *f_count = NULL;
+ hsize_t m_offset[1];
+ hsize_t m_count[1];
+ hsize_t msdims[1];
+ hsize_t size[3];
+ herr_t r;
+ herr_t check;
+ int check_pe;
+ int dim_flag;
+ int ndims;
+ int i;
+ char *subst_name = NULL;
+
+ if ( ev->path_length != 0 ) {
+ subst_name = partial_event_substitution(ev, name);
+ } else {
+ subst_name = strdup(name);
+ }
+
+ check_pe = check_path_existence(f->fh, subst_name);
+ if ( check_pe == 0 ) {
+ ERROR("No such event-based float field '%s'\n", subst_name);
+ return 1;
+ }
+
+ dh = H5Dopen2(f->fh, subst_name, H5P_DEFAULT);
+ type = H5Dget_type(dh);
+ class = H5Tget_class(type);
+
+ if ( class != H5T_FLOAT ) {
+ ERROR("Not a floating point value.\n");
+ H5Tclose(type);
+ H5Dclose(dh);
+ return 1;
+ }
+
+ /* Get the dimensionality. We have to cope with scalars expressed as
+ * arrays with all dimensions 1, as well as zero-d arrays. */
+ sh = H5Dget_space(dh);
+ ndims = H5Sget_simple_extent_ndims(sh);
+ if ( ndims > 3 ) {
+ H5Tclose(type);
+ H5Dclose(dh);
+ return 1;
+ }
+ H5Sget_simple_extent_dims(sh, size, NULL);
+
+ m_offset[0] = 0;
+ m_count[0] = 1;
+ msdims[0] = 1;
+ ms = H5Screate_simple(1,msdims,NULL);
+
+ /* Check that the size in all dimensions is 1 */
+ /* or that one of the dimensions has the same */
+ /* size as the hyperplane events */
+
+ dim_flag = 0;
+
+ for ( i=0; i<ndims; i++ ) {
+ if ( size[i] != 1 ) {
+ if ( i == 0 && size[i] > ev->dim_entries[0] ) {
+ dim_flag = 1;
+ } else {
+ H5Tclose(type);
+ H5Dclose(dh);
+ return 1;
+ }
+ }
+ }
+
+ if ( dim_flag == 0 ) {
+
+ r = H5Dread(dh, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, val);
+
+ if ( r < 0 ) {
+ ERROR("Couldn't read value.\n");
+ H5Tclose(type);
+ H5Dclose(dh);
+ return 1;
+ }
+
+ } else {
+
+ f_offset = malloc(ndims*sizeof(hsize_t));
+ f_count = malloc(ndims*sizeof(hsize_t));
+
+ for ( i=0; i<ndims; i++ ) {
+
+ if ( i == 0 ) {
+ f_offset[i] = ev->dim_entries[0];
+ f_count[i] = 1;
+ } else {
+ f_offset[i] = 0;
+ f_count[i] = 0;
+ }
+
+ }
+
+ check = H5Sselect_hyperslab(sh, H5S_SELECT_SET,
+ f_offset, NULL, f_count, NULL);
+ if ( check <0 ) {
+ ERROR("Error selecting dataspace for float value");
+ free(f_offset);
+ free(f_count);
+ return 1;
+ }
+
+ ms = H5Screate_simple(1,msdims,NULL);
+ check = H5Sselect_hyperslab(ms, H5S_SELECT_SET,
+ m_offset, NULL, m_count, NULL);
+ if ( check < 0 ) {
+ ERROR("Error selecting memory dataspace for float value");
+ free(f_offset);
+ free(f_count);
+ return 1;
+ }
+
+ r = H5Dread(dh, memtype, ms, sh, H5P_DEFAULT, val);
+ if ( r < 0 ) {
+ ERROR("Couldn't read value.\n");
+ H5Tclose(type);
+ H5Dclose(dh);
+ return 1;
+ }
+
+ }
+
+ free(subst_name);
+
+ return 0;
+}
+
+
+int hdfile_get_value(struct hdfile *f, const char *name,
+ struct event *ev, void *val, hid_t memtype)
+{
+ if ( ev == NULL ) {
+ return get_scalar_value(f, name, val, memtype);
+ } else {
+ return get_ev_based_value(f, name, ev, val, memtype);
+ }
+}
+
+
void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f,
struct event *ev, struct image *image)
{
@@ -1160,15 +1356,16 @@ void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f,
/* Explicit value given */
eV = beam->photon_energy;
- } else if ( ev != NULL ) {
-
- /* Value from HDF5 file, event-based structure */
- eV = get_ev_based_value(f, beam->photon_energy_from, ev);
-
} else {
- /* Value from HDF5 file, single-event structure */
- eV = get_value(f, beam->photon_energy_from);
+ int r;
+
+ r = hdfile_get_value(f, beam->photon_energy_from, ev, &eV,
+ H5T_NATIVE_DOUBLE);
+ if ( r ) {
+ ERROR("Failed to read '%s'\n",
+ beam->photon_energy_from);
+ }
}
@@ -1208,8 +1405,15 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element,
}
image->data = buf;
+ if ( image->det != NULL ) {
+ ERROR("WARNING: hdf5_read() called with geometry structure.\n");
+ }
+ image->det = simple_geometry(image);
+
if ( satcorr ) debodge_saturation(f, image);
+ unpack_panels(image, image->det);
+
if ( image->beam != NULL ) {
fill_in_beam_parameters(image->beam, f, NULL, image);
@@ -1588,255 +1792,6 @@ int hdfile_is_scalar(struct hdfile *f, const char *name, int verbose)
}
-
-
-static int get_f_value(struct hdfile *f, const char *name, double *val)
-{
- hid_t dh;
- hid_t type;
- hid_t class;
- herr_t r;
- double buf;
- int check;
-
- if ( !hdfile_is_scalar(f, name, 1) ) return 1;
-
- check = check_path_existence(f->fh, name);
- if ( check == 0 ) {
- ERROR("No such float field '%s'\n", name);
- return 1;
- }
-
- dh = H5Dopen2(f->fh, name, H5P_DEFAULT);
-
- type = H5Dget_type(dh);
- class = H5Tget_class(type);
-
- if ( class != H5T_FLOAT ) {
- ERROR("Not a floating point value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 1;
- }
-
- r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, &buf);
- if ( r < 0 ) {
- ERROR("Couldn't read value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 1;
- }
-
- *val = buf;
- return 0;
-}
-
-
-static int get_ev_based_f_value(struct hdfile *f, const char *name,
- struct event *ev, double *val)
-{
- hid_t dh;
- hid_t type;
- hid_t class;
- hid_t sh;
- hid_t ms;
- hsize_t *f_offset = NULL;
- hsize_t *f_count = NULL;
- hsize_t m_offset[1];
- hsize_t m_count[1];
- hsize_t msdims[1];
- hsize_t size[3];
- herr_t r;
- herr_t check;
- double buf;
- int check_pe;
- int dim_flag;
- int ndims;
- int i;
- char *subst_name = NULL;
-
- if ( ev->path_length != 0 ) {
- subst_name = partial_event_substitution(ev, name);
- } else {
- subst_name = strdup(name);
- }
-
- check_pe = check_path_existence(f->fh, subst_name);
- if ( check_pe == 0 ) {
- ERROR("No such event-based float field '%s'\n", subst_name);
- return 1;
- }
-
- dh = H5Dopen2(f->fh, name, H5P_DEFAULT);
- type = H5Dget_type(dh);
- class = H5Tget_class(type);
-
- if ( class != H5T_FLOAT ) {
- ERROR("Not a floating point value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 1;
- }
-
- /* Get the dimensionality. We have to cope with scalars expressed as
- * arrays with all dimensions 1, as well as zero-d arrays. */
- sh = H5Dget_space(dh);
- ndims = H5Sget_simple_extent_ndims(sh);
- if ( ndims > 3 ) {
- H5Tclose(type);
- H5Dclose(dh);
- return 1;
- }
- H5Sget_simple_extent_dims(sh, size, NULL);
-
- m_offset[0] = 0;
- m_count[0] = 1;
- msdims[0] = 1;
- ms = H5Screate_simple(1,msdims,NULL);
-
- /* Check that the size in all dimensions is 1 */
- /* or that one of the dimensions has the same */
- /* size as the hyperplane events */
-
- dim_flag = 0;
-
- for ( i=0; i<ndims; i++ ) {
- if ( size[i] != 1 ) {
- if ( i == 0 && size[i] > ev->dim_entries[0] ) {
- dim_flag = 1;
- } else {
- H5Tclose(type);
- H5Dclose(dh);
- return 1;
- }
- }
- }
-
- if ( dim_flag == 0 ) {
-
- r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, &buf);
-
- if ( r < 0 ) {
- ERROR("Couldn't read value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 1;
- }
-
- } else {
-
- f_offset = malloc(ndims*sizeof(hsize_t));
- f_count = malloc(ndims*sizeof(hsize_t));
-
- for ( i=0; i<ndims; i++ ) {
-
- if ( i == 0 ) {
- f_offset[i] = ev->dim_entries[0];
- f_count[i] = 1;
- } else {
- f_offset[i] = 0;
- f_count[i] = 0;
- }
-
- }
-
- check = H5Sselect_hyperslab(sh, H5S_SELECT_SET,
- f_offset, NULL, f_count, NULL);
- if ( check <0 ) {
- ERROR("Error selecting dataspace for float value");
- free(f_offset);
- free(f_count);
- return 1;
- }
-
- ms = H5Screate_simple(1,msdims,NULL);
- check = H5Sselect_hyperslab(ms, H5S_SELECT_SET,
- m_offset, NULL, m_count, NULL);
- if ( check < 0 ) {
- ERROR("Error selecting memory dataspace for float value");
- free(f_offset);
- free(f_count);
- return 1;
- }
-
- r = H5Dread(dh, H5T_NATIVE_DOUBLE, ms, sh,
- H5P_DEFAULT, &buf);
- if ( r < 0 ) {
- ERROR("Couldn't read value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 1;
- }
-
- }
-
- free(subst_name);
- *val = buf;
-
- return 0;
-}
-
-
-static int get_i_value(struct hdfile *f, const char *name, int *val)
-{
- hid_t dh;
- hid_t type;
- hid_t class;
- herr_t r;
- int buf;
- int check;
-
- if ( !hdfile_is_scalar(f, name, 1) ) return 1;
-
- check = check_path_existence(f->fh, name);
- if ( check == 0 ) {
- ERROR("No such integer field '%s'\n", name);
- return 1;
- }
-
- dh = H5Dopen2(f->fh, name, H5P_DEFAULT);
- type = H5Dget_type(dh);
- class = H5Tget_class(type);
-
- if ( class != H5T_INTEGER ) {
- ERROR("Not an integer value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 1;
- }
-
- r = H5Dread(dh, H5T_NATIVE_INT, H5S_ALL, H5S_ALL,
- H5P_DEFAULT, &buf);
- if ( r < 0 ) {
- ERROR("Couldn't read value.\n");
- H5Tclose(type);
- H5Dclose(dh);
- return 1;
- }
-
- *val = buf;
- return 0;
-}
-
-
-double get_value(struct hdfile *f, const char *name)
-{
- double val = 0.0;
- get_f_value(f, name, &val);
- return val;
-}
-
-double get_ev_based_value(struct hdfile *f, const char *name,
- struct event *ev)
-{
- double val = -1;
- get_ev_based_f_value(f, name, ev, &val);
- return val;
-}
-
-
struct copy_hdf5_field
{
char **fields;
@@ -1947,9 +1902,15 @@ char *hdfile_get_string_value(struct hdfile *f, const char *name,
hid_t class;
int buf_i;
double buf_f;
- char *tmp;
+ char *tmp = NULL, *subst_name = NULL;
- dh = H5Dopen2(f->fh, name, H5P_DEFAULT);
+ if (ev != NULL && ev->path_length != 0 ) {
+ subst_name = partial_event_substitution(ev, name);
+ } else {
+ subst_name = strdup(name);
+ }
+
+ dh = H5Dopen2(f->fh, subst_name, H5P_DEFAULT);
if ( dh < 0 ) return NULL;
type = H5Dget_type(dh);
class = H5Tget_class(type);
@@ -1957,7 +1918,6 @@ char *hdfile_get_string_value(struct hdfile *f, const char *name,
if ( class == H5T_STRING ) {
herr_t r;
- char *tmp;
hid_t sh;
size = H5Tget_size(type);
@@ -1966,51 +1926,54 @@ char *hdfile_get_string_value(struct hdfile *f, const char *name,
sh = H5Screate(H5S_SCALAR);
r = H5Dread(dh, type, sh, sh, H5P_DEFAULT, tmp);
- if ( r < 0 ) goto fail;
-
- /* Two possibilities:
- * String is already zero-terminated
- * String is not terminated.
- * Make sure things are done properly... */
- tmp[size] = '\0';
- chomp(tmp);
+ if ( r < 0 ) {
+ free(tmp);
+ tmp = NULL;
+ } else {
- return tmp;
+ /* Two possibilities:
+ * String is already zero-terminated
+ * String is not terminated.
+ * Make sure things are done properly... */
+ tmp[size] = '\0';
+ chomp(tmp);
+ }
+ } else {
- }
+ int r;
- switch ( class ) {
+ switch ( class ) {
- case H5T_FLOAT :
- if ( ev != NULL ) {
- if ( get_ev_based_f_value(f, name, ev, &buf_f) ) goto fail;
- } else {
- if ( get_f_value(f, name, &buf_f) ) goto fail;
+ case H5T_FLOAT :
+ r = hdfile_get_value(f, subst_name, ev, &buf_f,
+ H5T_NATIVE_DOUBLE);
+ if ( r == 0 ) {
+ tmp = malloc(256);
+ snprintf(tmp, 255, "%f", buf_f);
+ }
+ break;
+
+ case H5T_INTEGER :
+ r = hdfile_get_value(f, subst_name, ev, &buf_i,
+ H5T_NATIVE_INT);
+ if ( r == 0 ) {
+ tmp = malloc(256);
+ snprintf(tmp, 255, "%d", buf_i);
+ }
+ break;
}
- tmp = malloc(256);
- snprintf(tmp, 255, "%f", buf_f);
- return tmp;
-
- case H5T_INTEGER :
- if ( get_i_value(f, name, &buf_i) ) goto fail;
- tmp = malloc(256);
- snprintf(tmp, 255, "%d", buf_i);
- return tmp;
-
- default :
- goto fail;
}
-fail:
H5Tclose(type);
H5Dclose(dh);
- return NULL;
+ free(subst_name);
+ return tmp;
}
char **hdfile_read_group(struct hdfile *f, int *n, const char *parent,
- int **p_is_group, int **p_is_image)
+ int **p_is_group, int **p_is_image)
{
hid_t gh;
hsize_t num;
diff --git a/libcrystfel/src/hdf5-file.h b/libcrystfel/src/hdf5-file.h
index 4d7b223c..8c89eb93 100644
--- a/libcrystfel/src/hdf5-file.h
+++ b/libcrystfel/src/hdf5-file.h
@@ -3,11 +3,11 @@
*
* Read/write HDF5 data files
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2009-2012 Thomas White <taw@physics.org>
+ * 2009-2015 Thomas White <taw@physics.org>
* 2014 Valerio Mariani
*
@@ -76,15 +76,10 @@ extern char **hdfile_read_group(struct hdfile *f, int *n, const char *parent,
extern int hdfile_set_first_image(struct hdfile *f, const char *group);
extern void hdfile_close(struct hdfile *f);
-extern int hdfile_is_scalar(struct hdfile *f, const char *name, int verbose);
-char *hdfile_get_string_value(struct hdfile *f, const char *name,
- struct event* ev);
-extern int get_peaks(struct image *image, struct hdfile *f, const char *p,
- int cxi_format, struct filename_plus_event *fpe);
-extern double get_value(struct hdfile *f, const char *name);
+extern int get_peaks(struct image *image, struct hdfile *f, const char *p);
-extern double get_ev_based_value(struct hdfile *f, const char *name,
- struct event *ev);
+extern int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
+ struct filename_plus_event *fpe);
extern struct copy_hdf5_field *new_copy_hdf5_field_list(void);
extern void free_copy_hdf5_field_list(struct copy_hdf5_field *f);
@@ -97,6 +92,12 @@ extern void add_copy_hdf5_field(struct copy_hdf5_field *copyme,
extern struct event_list *fill_event_list(struct hdfile* hdfile,
struct detector* det);
+extern int hdfile_get_value(struct hdfile *f, const char *name,
+ struct event *ev, void *val, hid_t memtype);
+extern int hdfile_is_scalar(struct hdfile *f, const char *name, int verbose);
+extern char *hdfile_get_string_value(struct hdfile *f, const char *name,
+ struct event *ev);
+
#ifdef __cplusplus
}
#endif
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 97187655..6b36af3e 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -56,6 +56,19 @@
#include "grainspotter.h"
+static int debug_index(struct image *image)
+{
+ Crystal *cr = crystal_new();
+ UnitCell *cell = cell_new();
+ cell_set_reciprocal(cell, +0.0000e9, +0.0000e9, +0.0000e9,
+ +0.0000e9, +0.0000e9, +0.0000e9,
+ +0.0000e9, +0.0000e9, +0.0000e9);
+ crystal_set_cell(cr, cell);
+ image_add_crystal(image, cr);
+ return 1;
+}
+
+
IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
struct detector *det, float *ltl)
{
@@ -101,6 +114,10 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
det, ltl);
break;
+ case INDEXING_DEBUG :
+ iprivs[n] = (IndexingPrivate *)strdup("Hello!");
+ break;
+
default :
ERROR("Don't know how to prepare indexing method %i\n",
indm[n]);
@@ -175,6 +192,10 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs)
grainspotter_cleanup(privs[n]);
break;
+ case INDEXING_DEBUG :
+ free(privs[n]);
+ break;
+
default :
ERROR("Don't know how to clean up indexing method %i\n",
indms[n]);
@@ -244,6 +265,9 @@ static int try_indexer(struct image *image, IndexingMethod indm,
return grainspotter_index(image, ipriv);
break;
+ case INDEXING_DEBUG :
+ return debug_index(image);
+
default :
ERROR("Unrecognised indexing method: %i\n", indm);
break;
@@ -381,8 +405,12 @@ char *indexer_str(IndexingMethod indm)
strcpy(str, "simulation");
break;
+ case INDEXING_DEBUG :
+ strcpy(str, "debug");
+ break;
+
default :
- ERROR("Unrecognised indexing method %i\n",
+ ERROR("No test description for indexing method %i\n",
indm & INDEXING_METHOD_MASK);
strcpy(str, "(unknown)");
break;
@@ -457,6 +485,10 @@ IndexingMethod *build_indexer_list(const char *str)
list[++nmeth] = INDEXING_SIMULATION;
return list;
+ } else if ( strcmp(methods[i], "debug") == 0) {
+ list[++nmeth] = INDEXING_DEBUG;
+ return list;
+
} else if ( strcmp(methods[i], "raw") == 0) {
list[nmeth] = set_raw(list[nmeth]);
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index 709b7507..ff1cd5df 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -3,13 +3,13 @@
*
* Perform indexing (somehow)
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2010-2014 Thomas White <taw@physics.org>
+ * 2010-2015 Thomas White <taw@physics.org>
* 2010 Richard Kirian
* 2012 Lorenzo Galli
*
@@ -72,6 +72,7 @@
* @INDEXING_GRAINSPOTTER: Invoke GrainSpotter
* @INDEXING_XDS: Invoke XDS
* @INDEXING_SIMULATION: Dummy value
+ * @INDEXING_DEBUG: Results injector for debugging
* @INDEXING_CHECK_CELL_COMBINATIONS: Check linear combinations of unit cell
* axes for agreement with given cell.
* @INDEXING_CHECK_CELL_AXES: Check unit cell axes for agreement with given
@@ -98,7 +99,8 @@ typedef enum {
INDEXING_GRAINSPOTTER = 4,
INDEXING_XDS = 5,
INDEXING_SIMULATION = 6,
- INDEXING_ASDF = 7,
+ INDEXING_DEBUG = 7,
+ INDEXING_ASDF = 8,
/* Bits at the top of the IndexingMethod are flags which modify the
* behaviour of the indexer. */
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index 4a743c63..16c27bbe 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -54,94 +54,6 @@
#include "integration.h"
-static void check_eigen(gsl_vector *e_val)
-{
- int i;
- double vmax, vmin;
- const int n = e_val->size;
- const double max_condition = 1e6;
- const int verbose = 0;
-
- if ( verbose ) STATUS("Eigenvalues:\n");
- vmin = +INFINITY;
- vmax = 0.0;
- for ( i=0; i<n; i++ ) {
- double val = gsl_vector_get(e_val, i);
- if ( verbose ) STATUS("%i: %e\n", i, val);
- if ( val > vmax ) vmax = val;
- if ( val < vmin ) vmin = val;
- }
-
- for ( i=0; i<n; i++ ) {
- double val = gsl_vector_get(e_val, i);
- if ( val < vmax/max_condition ) {
- gsl_vector_set(e_val, i, 0.0);
- }
- }
-
- vmin = +INFINITY;
- vmax = 0.0;
- for ( i=0; i<n; i++ ) {
- double val = gsl_vector_get(e_val, i);
- if ( val == 0.0 ) continue;
- if ( val > vmax ) vmax = val;
- if ( val < vmin ) vmin = val;
- }
- if ( verbose ) {
- STATUS("Condition number: %e / %e = %5.2f\n",
- vmax, vmin, vmax/vmin);
- }
-}
-
-
-static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *Mp)
-{
- gsl_matrix *s_vec;
- gsl_vector *s_val;
- int err, n;
- gsl_vector *shifts;
- gsl_matrix *M;
-
- n = v->size;
- if ( v->size != Mp->size1 ) return NULL;
- if ( v->size != Mp->size2 ) return NULL;
-
- M = gsl_matrix_alloc(n, n);
- if ( M == NULL ) return NULL;
- gsl_matrix_memcpy(M, Mp);
-
- s_val = gsl_vector_calloc(n);
- s_vec = gsl_matrix_calloc(n, n);
-
- err = gsl_linalg_SV_decomp_jacobi(M, s_vec, s_val);
- if ( err ) {
- ERROR("SVD failed: %s\n", gsl_strerror(err));
- gsl_matrix_free(s_vec);
- gsl_vector_free(s_val);
- return NULL;
- }
- /* "M" is now "U" */
-
- check_eigen(s_val);
-
- shifts = gsl_vector_calloc(n);
- err = gsl_linalg_SV_solve(M, s_vec, s_val, v, shifts);
- if ( err ) {
- ERROR("Matrix solution failed: %s\n", gsl_strerror(err));
- gsl_matrix_free(s_vec);
- gsl_vector_free(s_val);
- gsl_vector_free(shifts);
- return NULL;
- }
-
- gsl_matrix_free(s_vec);
- gsl_vector_free(s_val);
- gsl_matrix_free(M);
-
- return shifts;
-}
-
-
enum boxmask_val
{
BM_IG, /* "Soft" ignore */
@@ -340,14 +252,14 @@ static void show_reference_profile(struct intcontext *ic, int i)
static void show_peak_box(struct intcontext *ic, struct peak_box *bx,
- int results_pipe)
+ pthread_mutex_t *term_lock)
{
#ifdef HAVE_CURSES_COLOR
int q;
signed int h, k, l;
double fs, ss;
- if ( results_pipe != 0 ) write(results_pipe, "SUSPEND\n", 8);
+ if ( term_lock != NULL ) pthread_mutex_lock(term_lock);
initscr();
clear();
@@ -411,7 +323,10 @@ static void show_peak_box(struct intcontext *ic, struct peak_box *bx,
getch();
endwin();
- if ( results_pipe != 0 ) write(results_pipe, "RELEASE\n", 8);
+ if ( term_lock != NULL ) pthread_mutex_unlock(term_lock);
+#else
+ STATUS("Not showing peak box because CrystFEL was compiled without "
+ "ncurses.\n");
#endif
}
@@ -442,8 +357,8 @@ static void fit_bg(struct intcontext *ic, struct peak_box *bx)
}
}
- /* SVD is massive overkill here */
- ans = solve_svd(v, bx->bgm);
+ /* FIXME: SVD is massive overkill here */
+ ans = solve_svd(v, bx->bgm, NULL, 0);
gsl_vector_free(v);
bx->a = gsl_vector_get(ans, 0);
@@ -500,7 +415,7 @@ static int init_intcontext(struct intcontext *ic)
}
/* How many reference profiles? */
- ic->n_reference_profiles = ic->image->det->n_panels;
+ ic->n_reference_profiles = 1;
ic->reference_profiles = calloc(ic->n_reference_profiles,
sizeof(double *));
if ( ic->reference_profiles == NULL ) return 1;
@@ -1255,7 +1170,7 @@ static int get_int_diag(struct intcontext *ic, Reflection *refl)
static void integrate_prof2d_once(struct intcontext *ic, struct peak_box *bx,
- int results_pipe)
+ pthread_mutex_t *term_lock)
{
bx->intensity = fit_intensity(ic, bx);
bx->sigma = calc_sigma(ic, bx);
@@ -1278,15 +1193,16 @@ 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++;
set_redundancy(bx->refl, 0);
}
- if ( get_int_diag(ic, bx->refl) ) show_peak_box(ic, bx,
- results_pipe);
+ if ( get_int_diag(ic, bx->refl) ) {
+ show_peak_box(ic, bx, term_lock);
+ }
} else {
@@ -1379,7 +1295,7 @@ static void integrate_prof2d(IntegrationMethod meth,
Crystal *cr, struct image *image, IntDiag int_diag,
signed int idh, signed int idk, signed int idl,
double ir_inn, double ir_mid, double ir_out,
- int results_pipe, int **masks)
+ pthread_mutex_t *term_lock, int **masks)
{
RefList *list;
UnitCell *cell;
@@ -1423,7 +1339,7 @@ static void integrate_prof2d(IntegrationMethod meth,
for ( i=0; i<ic.n_boxes; i++ ) {
struct peak_box *bx;
bx = &ic.boxes[i];
- integrate_prof2d_once(&ic, bx, results_pipe);
+ integrate_prof2d_once(&ic, bx, term_lock);
}
//refine_rigid_groups(&ic);
@@ -1436,7 +1352,7 @@ static void integrate_prof2d(IntegrationMethod meth,
static void integrate_rings_once(Reflection *refl, struct image *image,
struct intcontext *ic, UnitCell *cell,
- int results_pipe)
+ pthread_mutex_t *term_lock)
{
double pfs, pss;
struct peak_box *bx;
@@ -1482,6 +1398,7 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
r = check_box(ic, bx, &saturated);
if ( !r ) {
fit_bg(ic, bx);
+ if ( !bg_ok(bx) ) r = 1;
}
bx->offs_fs = 0.0;
bx->offs_ss = 0.0;
@@ -1532,9 +1449,9 @@ 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);
+ if ( get_int_diag(ic, refl) ) show_peak_box(ic, bx, term_lock);
if ( intensity < -5.0*sigma ) {
ic->n_implausible++;
@@ -1634,7 +1551,7 @@ static void integrate_rings(IntegrationMethod meth,
Crystal *cr, struct image *image, IntDiag int_diag,
signed int idh, signed int idk, signed int idl,
double ir_inn, double ir_mid, double ir_out,
- int results_pipe, int **masks)
+ pthread_mutex_t *term_lock, int **masks)
{
RefList *list;
Reflection *refl;
@@ -1670,7 +1587,7 @@ static void integrate_rings(IntegrationMethod meth,
refl != NULL;
refl = next_refl(refl, iter) )
{
- integrate_rings_once(refl, image, &ic, cell, results_pipe);
+ integrate_rings_once(refl, image, &ic, cell, term_lock);
}
//refine_rigid_groups(&ic);
@@ -1687,7 +1604,7 @@ void integrate_all_4(struct image *image, IntegrationMethod meth,
double ir_inn, double ir_mid, double ir_out,
IntDiag int_diag,
signed int idh, signed int idk, signed int idl,
- int results_pipe)
+ pthread_mutex_t *term_lock)
{
int i;
int *masks[image->det->n_panels];
@@ -1727,14 +1644,14 @@ void integrate_all_4(struct image *image, IntegrationMethod meth,
integrate_rings(meth, cr, image,
int_diag, idh, idk, idl,
ir_inn, ir_mid, ir_out,
- results_pipe, masks);
+ term_lock, masks);
break;
case INTEGRATION_PROF2D :
integrate_prof2d(meth, cr, image,
int_diag, idh, idk, idl,
ir_inn, ir_mid, ir_out,
- results_pipe, masks);
+ term_lock, masks);
break;
default :
diff --git a/libcrystfel/src/integration.h b/libcrystfel/src/integration.h
index ca04157a..a012fc14 100644
--- a/libcrystfel/src/integration.h
+++ b/libcrystfel/src/integration.h
@@ -3,11 +3,11 @@
*
* Integration of intensities
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2014 Thomas White <taw@physics.org>
+ * 2010-2015 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -33,7 +33,6 @@
#include <config.h>
#endif
-
#include "geometry.h"
/**
@@ -127,7 +126,7 @@ extern void integrate_all_4(struct image *image, IntegrationMethod meth,
double ir_inn, double ir_mid, double ir_out,
IntDiag int_diag,
signed int idh, signed int idk, signed int idl,
- int results_pipe);
+ pthread_mutex_t *term_lock);
#ifdef __cplusplus
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index 9900c232..ba724419 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -3,11 +3,11 @@
*
* Peak search and other image analysis
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2014 Thomas White <taw@physics.org>
+ * 2010-2015 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
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..ccc421f3 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;
@@ -105,6 +102,7 @@ struct _refldata {
/* User-specified temporary values */
double temp1;
double temp2;
+ int flag;
};
@@ -316,29 +314,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;
}
@@ -547,6 +547,22 @@ double get_temp2(const Reflection *refl)
}
+/**
+ * get_flag:
+ * @refl: A %Reflection
+ *
+ * The integer flag value can be used according to the needs of the calling
+ * program.
+ *
+ * Returns: the flag for this reflection.
+ *
+ **/
+int get_flag(const Reflection *refl)
+{
+ return refl->data.flag;
+}
+
+
/********************************** Setters ***********************************/
/**
@@ -570,20 +586,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
@@ -761,6 +789,21 @@ void set_temp2(Reflection *refl, double temp)
}
+/**
+ * set_flag
+ * @refl: A %Reflection
+ * @flag: New flag value
+ *
+ * @flag is an integer value which can be used according to the needs of the
+ * calling program.
+ *
+ **/
+void set_flag(Reflection *refl, int flag)
+{
+ refl->data.flag = flag;
+}
+
+
/********************************* Insertion **********************************/
static Reflection *rotate_once(Reflection *refl, int dir)
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index 85a87c54..dac313a4 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,
@@ -103,11 +103,12 @@ extern double get_esd_intensity(const Reflection *refl);
extern double get_phase(const Reflection *refl, int *have_phase);
extern double get_peak(const Reflection *refl);
extern double get_mean_bg(const Reflection *refl);
+extern int get_flag(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);
@@ -121,6 +122,7 @@ extern void set_peak(Reflection *refl, double peak);
extern void set_mean_bg(Reflection *refl, double mean_bg);
extern void set_symmetric_indices(Reflection *refl,
signed int hs, signed int ks, signed int ls);
+extern void set_flag(Reflection *refl, int flag);
/* Insertion */
extern Reflection *add_refl(RefList *list,
diff --git a/libcrystfel/src/render.c b/libcrystfel/src/render.c
index 17318782..2dcb7b93 100644
--- a/libcrystfel/src/render.c
+++ b/libcrystfel/src/render.c
@@ -104,6 +104,84 @@ static void render_rgb(double val, double max,
}
+static void render_geoptimiser(double val, double max,
+ double *rp, double *gp, double *bp)
+{
+ double r;
+ double p;
+
+ r = val/max;
+
+ if ( val < 0.0 ) {
+ *rp = 0.0;
+ *gp = 0.0;
+ *bp = 0.0;
+ return;
+ }
+
+ if ( r >= 0.0 && r < 0.059 ) {
+ p = (r-0.0)/(0.059-0.0);
+ *rp = 0.0;
+ *gp = 0.0;
+ *bp = ((91.0/256.0)-0.0)*p;
+ return;
+ }
+
+ if ( r >= 0.059 && r < 0.220 ) {
+ p = (r-0.059)/(0.220-0.059);
+ *rp = ((122.0/256.0)-0.0)*p;
+ *gp = 0.0;
+ *bp = ((227.0/256.0)-(91.0/256.0))*p+(91.0/256.0);
+ return;
+ }
+
+ if ( r >= 0.220 && r < 0.376 ) {
+ p = (r-0.220)/(0.376-0.220);
+ *rp = ((195.0/256.0)-(122.0/256.0))*p+(122.0/256.0);
+ *gp = 0.0;
+ *bp = ((93.0/256.0)-(227.0/256.0))*p+(227.0/256.0);
+ return;
+ }
+
+ if ( r >= 0.376 && r < 0.498 ) {
+ p = (r-0.376)/(0.498-0.376);
+ *rp = ((238.0/256.0)-(195.0/256.0))*p+(195.0/256.0);
+ *gp = ((76.0/256.0)-0.0)*p;
+ *bp = (0.0-(93.0/256.0))*p+(93.0/256.0);
+ return;
+ }
+
+ if ( r >= 0.498 && r < 0.564 ) {
+ p = (r-0.498)/(0.564-0.498);
+ *rp = (1.0-(238.0/256.0))*p+(238.0/256.0);
+ *gp = ((117.0/256.0)-(76.0/256.0))*p+(76.0/256.0);
+ *bp = 0.0;
+ return;
+ }
+
+ if ( r >= 0.564 && r < 0.815 ) {
+ p = (r-0.564)/(0.815-0.564);
+ *rp = 1.0;
+ *gp = ((234.0/256.0)-(117.0/256.0))*p+(117.0/256.0);
+ *bp = 0.0;
+ return;
+ }
+
+ if ( r >= 0.815 && r < 1.0 ) {
+ p = (r-0.815)/(1.0-0.815);
+ *rp = 1.0;
+ *gp = (1.0-(234.0/256.0))*p+(234.0/256.0);
+ *bp = (1.0-0.0)*p;
+ return;
+ }
+
+ if ( r >= 1.0 ) {
+ *rp = 1.0; *gp = 1.0; *bp = 1.0;
+ return;
+ }
+}
+
+
static void render_ratio(double val, double max,
double *rp, double *gp, double *bp)
{
@@ -164,5 +242,9 @@ void render_scale(double val, double max, int scale,
case SCALE_RATIO :
render_ratio(val, max, rp, gp, bp);
break;
+
+ case SCALE_GEOPTIMISER :
+ render_geoptimiser(val, max, rp, gp, bp);
+ break;
}
}
diff --git a/libcrystfel/src/render.h b/libcrystfel/src/render.h
index a4da922d..c2e5965a 100644
--- a/libcrystfel/src/render.h
+++ b/libcrystfel/src/render.h
@@ -38,7 +38,8 @@ enum {
SCALE_COLOUR,
SCALE_MONO,
SCALE_INVMONO,
- SCALE_RATIO
+ SCALE_RATIO,
+ SCALE_GEOPTIMISER
};
#ifdef __cplusplus
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 118cd782..808149e3 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -3,12 +3,12 @@
*
* Stream tools
*
- * 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>
* 2014 Valerio Mariani
* 2011 Richard Kirian
* 2011 Andrew Aquila
@@ -319,7 +319,8 @@ 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_panel(refl, p);
}
set_esd_intensity(refl, sigma);
set_peak(refl, pk);
@@ -383,11 +384,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 +449,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);
}
@@ -708,6 +709,10 @@ static int write_crystal(Stream *st, Crystal *cr, int include_reflections)
rad = crystal_get_profile_radius(cr);
fprintf(st->fh, "profile_radius = %.5f nm^-1\n", rad/1e9);
+ if ( crystal_get_notes(cr) != NULL ) {
+ fprintf(st->fh, "%s\n", crystal_get_notes(cr));
+ }
+
reflist = crystal_get_reflections(cr);
if ( reflist != NULL ) {
@@ -838,7 +843,10 @@ int write_chunk(Stream *st, struct image *i, struct hdfile *hdfile,
}
for ( j=0; j<i->n_crystals; j++ ) {
- ret = write_crystal(st, i->crystals[j], include_reflections);
+ if ( crystal_get_user_flag(i->crystals[j]) == 0 ) {
+ ret = write_crystal(st, i->crystals[j],
+ include_reflections);
+ }
}
fprintf(st->fh, CHUNK_END_MARKER"\n");
@@ -1260,7 +1268,12 @@ Stream *open_stream_for_read(const char *filename)
st = malloc(sizeof(struct _stream));
if ( st == NULL ) return NULL;
- st->fh = fopen(filename, "r");
+ if ( strcmp(filename, "-") == 0 ) {
+ st->fh = stdin;
+ } else {
+ st->fh = fopen(filename, "r");
+ }
+
if ( st->fh == NULL ) {
free(st);
return NULL;
@@ -1499,7 +1512,7 @@ void write_geometry_file(Stream *st, const char *geom_filename) {
do {
rval = fgets(line, 1023, geom_fh);
- fputs(line, st->fh);
+ if ( rval != NULL ) fputs(line, st->fh);
} while ( rval != NULL );
fclose(geom_fh);
@@ -1526,27 +1539,3 @@ int rewind_stream(Stream *st)
return fseek(st->fh, 0, SEEK_SET);
}
-
-
-double extract_f_from_stuff(const char *field_name,
- struct stuff_from_stream* stuff)
-{
- int i;
-
- char field_name_plus_equal[256];
- sprintf(field_name_plus_equal, "hdf5%s = ", field_name);
-
-
-
- for ( i=0; i<stuff->n_fields; i++ ) {
-
- if ( strncmp(stuff->fields[i], field_name_plus_equal,
- strlen(field_name_plus_equal)) == 0 ) {
- return atoi(stuff->fields[i]+
- strlen(field_name_plus_equal));
- }
- }
-
- ERROR("Failed to recovery camera length from stream file\n");
- return -1;
-}
diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h
index afa9acda..a8cf4639 100644
--- a/libcrystfel/src/stream.h
+++ b/libcrystfel/src/stream.h
@@ -109,9 +109,6 @@ extern void write_geometry_file(Stream *st, const char *geom_filename);
extern int rewind_stream(Stream *st);
extern int is_stream(const char *filename);
-extern double extract_f_from_stuff(const char *field_name,
- struct stuff_from_stream* stuff);
-
#ifdef __cplusplus
}
#endif
diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c
index f8f93046..af0abedc 100644
--- a/libcrystfel/src/utils.c
+++ b/libcrystfel/src/utils.c
@@ -40,6 +40,8 @@
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_eigen.h>
#include "utils.h"
#include "image.h"
@@ -108,6 +110,161 @@ void show_matrix(gsl_matrix *M)
}
+static int check_eigen(gsl_vector *e_val, int verbose)
+{
+ int i;
+ double vmax, vmin;
+ const int n = e_val->size;
+ const double max_condition = 1e6;
+ int n_filt = 0;
+
+ if ( verbose ) STATUS("Eigenvalues:\n");
+ vmin = +INFINITY;
+ vmax = 0.0;
+ for ( i=0; i<n; i++ ) {
+ double val = gsl_vector_get(e_val, i);
+ if ( verbose ) STATUS("%i: %e\n", i, val);
+ if ( val > vmax ) vmax = val;
+ if ( val < vmin ) vmin = val;
+ }
+
+ for ( i=0; i<n; i++ ) {
+ double val = gsl_vector_get(e_val, i);
+ if ( val < vmax/max_condition ) {
+ gsl_vector_set(e_val, i, 0.0);
+ n_filt++;
+ }
+ }
+
+ vmin = +INFINITY;
+ vmax = 0.0;
+ for ( i=0; i<n; i++ ) {
+ double val = gsl_vector_get(e_val, i);
+ if ( val == 0.0 ) continue;
+ if ( val > vmax ) vmax = val;
+ if ( val < vmin ) vmin = val;
+ }
+ if ( verbose ) {
+ STATUS("Condition number: %e / %e = %5.2f\n",
+ vmax, vmin, vmax/vmin);
+ STATUS("%i out of %i eigenvalues filtered.\n", n_filt, n);
+ }
+
+ return n_filt;
+}
+
+
+/**
+ * solve_svd:
+ * v: a gsl_vector
+ * M: a gsl_matrix
+ * n_filt: pointer to store the number of filtered eigenvalues
+ * verbose: flag for verbosity on the terminal
+ *
+ * Solves the matrix equation M.x = v, returning x.
+ * Performs rescaling and eigenvalue filtering.
+ **/
+gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *pn_filt, int verbose)
+{
+ gsl_matrix *s_vec;
+ gsl_vector *s_val;
+ int err, n;
+ gsl_vector *shifts;
+ gsl_vector *SB;
+ gsl_vector *SinvX;
+ gsl_matrix *S; /* rescaling matrix due to Bricogne */
+ gsl_matrix *AS;
+ gsl_matrix *SAS;
+ int i;
+ int n_filt;
+ gsl_matrix *SAS_copy;
+
+ n = v->size;
+ if ( v->size != M->size1 ) return NULL;
+ if ( v->size != M->size2 ) return NULL;
+
+ /* Calculate the rescaling matrix S */
+ S = gsl_matrix_calloc(n, n);
+ for ( i=0; i<n; i++ ) {
+ double sii = pow(gsl_matrix_get(M, i, i), -0.5);
+ gsl_matrix_set(S, i, i, sii);
+ }
+
+ /* Calculate the matrix SAS, which we will be (not) inverting */
+ AS = gsl_matrix_calloc(n, n);
+ SAS = gsl_matrix_calloc(n, n);
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, M, S, 0.0, AS);
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, S, AS, 0.0, SAS);
+ gsl_matrix_free(AS);
+
+ /* Calculate the vector SB, which is the RHS of the equation */
+ SB = gsl_vector_calloc(n);
+ gsl_blas_dgemv(CblasNoTrans, 1.0, S, v, 0.0, SB);
+
+ if ( verbose ) {
+ STATUS("The equation after rescaling:\n");
+ show_matrix_eqn(SAS, SB);
+ }
+
+ SAS_copy = gsl_matrix_alloc(SAS->size1, SAS->size2);
+ gsl_matrix_memcpy(SAS_copy, SAS);
+
+ for ( i=0; i<n; i++ ) {
+ int j;
+ if ( isnan(gsl_vector_get(SB, i)) ) gsl_vector_set(SB, i, 0.0);
+ for ( j=0; j<n; j++ ) {
+ if ( isnan(gsl_matrix_get(SAS, i, j)) ) {
+ gsl_matrix_set(SAS, i, j, 0.0);
+ }
+ }
+ }
+
+ /* Do the SVD */
+ s_val = gsl_vector_calloc(n);
+ s_vec = gsl_matrix_calloc(n, n);
+ err = gsl_linalg_SV_decomp_jacobi(SAS, s_vec, s_val);
+ if ( err ) {
+ if ( verbose ) ERROR("SVD failed: %s\n", gsl_strerror(err));
+ gsl_matrix_free(s_vec);
+ gsl_vector_free(s_val);
+ gsl_matrix_free(SAS);
+ gsl_matrix_free(S);
+ return NULL;
+ }
+ /* "SAS" is now "U" */
+
+ /* Filter the eigenvalues */
+ n_filt = check_eigen(s_val, verbose);
+ if ( pn_filt != NULL ) *pn_filt = n_filt;
+
+ gsl_matrix_free(SAS_copy);
+
+ /* Solve the equation SAS.SinvX = SB */
+ SinvX = gsl_vector_calloc(n);
+ err = gsl_linalg_SV_solve(SAS, s_vec, s_val, SB, SinvX);
+ gsl_vector_free(SB);
+ gsl_matrix_free(SAS);
+ gsl_matrix_free(s_vec);
+ gsl_vector_free(s_val);
+
+ if ( err ) {
+ ERROR("Matrix solution failed: %s\n", gsl_strerror(err));
+ gsl_matrix_free(S);
+ gsl_vector_free(SinvX);
+ return NULL;
+ }
+
+ /* Calculate S.SinvX to get X, the shifts */
+ shifts = gsl_vector_calloc(n);
+ gsl_blas_dgemv(CblasNoTrans, 1.0, S, SinvX, 0.0, shifts);
+
+ gsl_matrix_free(S);
+ gsl_vector_free(SinvX);
+
+ return shifts;
+}
+
+
size_t notrail(char *s)
{
ssize_t i;
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index 2f2009d6..4955f875 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -103,6 +103,8 @@ extern struct rvec quat_rot(struct rvec q, struct quaternion z);
extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v);
extern void show_matrix(gsl_matrix *M);
+extern gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt,
+ int verbose);
extern size_t notrail(char *s);
extern void chomp(char *s);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 22:24:55 +0000