aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/asdf.c16
-rw-r--r--libcrystfel/src/asdf.h24
-rw-r--r--libcrystfel/src/cell-utils.c6
-rw-r--r--libcrystfel/src/cell.c10
-rw-r--r--libcrystfel/src/detector.c78
-rw-r--r--libcrystfel/src/detector.h49
-rw-r--r--libcrystfel/src/dirax.c12
-rw-r--r--libcrystfel/src/dirax.h13
-rw-r--r--libcrystfel/src/events.c57
-rw-r--r--libcrystfel/src/events.h1
-rw-r--r--libcrystfel/src/felix.c9
-rw-r--r--libcrystfel/src/felix.h16
-rw-r--r--libcrystfel/src/hdf5-file.c217
-rw-r--r--libcrystfel/src/hdf5-file.h15
-rw-r--r--libcrystfel/src/image.c698
-rw-r--r--libcrystfel/src/image.h134
-rw-r--r--libcrystfel/src/index.c416
-rw-r--r--libcrystfel/src/index.h27
-rw-r--r--libcrystfel/src/mosflm.c8
-rw-r--r--libcrystfel/src/mosflm.h12
-rw-r--r--libcrystfel/src/peakfinder8.c1195
-rw-r--r--libcrystfel/src/peakfinder8.h50
-rw-r--r--libcrystfel/src/peaks.c28
-rw-r--r--libcrystfel/src/peaks.h14
-rw-r--r--libcrystfel/src/reflist-utils.c2
-rw-r--r--libcrystfel/src/statistics.c1
-rw-r--r--libcrystfel/src/stream.c14
-rw-r--r--libcrystfel/src/stream.h12
-rw-r--r--libcrystfel/src/utils.c276
-rw-r--r--libcrystfel/src/utils.h86
-rw-r--r--libcrystfel/src/xds.c14
-rw-r--r--libcrystfel/src/xds.h16
32 files changed, 2926 insertions, 600 deletions
diff --git a/libcrystfel/src/asdf.c b/libcrystfel/src/asdf.c
index 49917ad9..0c43fe7a 100644
--- a/libcrystfel/src/asdf.c
+++ b/libcrystfel/src/asdf.c
@@ -4,12 +4,12 @@
* Alexandra's Superior Direction Finder, or
* Algorithm Similar to DirAx, FFT-based
*
- * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
* 2014-2015 Alexandra Tolstikova <alexandra.tolstikova@desy.de>
- * 2015 Thomas White <taw@physics.org>
+ * 2015,2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -28,8 +28,6 @@
*
*/
-#define _ISOC99_SOURCE
-
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
@@ -1103,7 +1101,7 @@ static int index_refls(gsl_vector **reflections, int N_reflections,
}
-int run_asdf(struct image *image, IndexingPrivate *ipriv)
+int run_asdf(struct image *image, void *ipriv)
{
int i, j;
@@ -1204,8 +1202,8 @@ int run_asdf(struct image *image, IndexingPrivate *ipriv)
}
-IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell,
- struct detector *det, float *ltl)
+void *asdf_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct detector *det, float *ltl)
{
struct asdf_private *dp;
int need_cell = 0;
@@ -1234,11 +1232,11 @@ IndexingPrivate *asdf_prepare(IndexingMethod *indm, UnitCell *cell,
dp->fftw = fftw_vars_new();
- return (IndexingPrivate *)dp;
+ return (void *)dp;
}
-void asdf_cleanup(IndexingPrivate *pp)
+void asdf_cleanup(void *pp)
{
struct asdf_private *p;
p = (struct asdf_private *)pp;
diff --git a/libcrystfel/src/asdf.h b/libcrystfel/src/asdf.h
index 1e492a6f..f130d63d 100644
--- a/libcrystfel/src/asdf.h
+++ b/libcrystfel/src/asdf.h
@@ -4,12 +4,12 @@
* Alexandra's Superior Direction Finder, or
* Algorithm Similar to DirAx, FFT-based
*
- * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
* 2014-2015 Alexandra Tolstikova <alexandra.tolstikova@desy.de>
- * 2015 Thomas White <taw@physics.org>
+ * 2015,2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -43,33 +43,33 @@ extern "C" {
#ifdef HAVE_FFTW
-extern int run_asdf(struct image *image, IndexingPrivate *ipriv);
+extern int run_asdf(struct image *image, void *ipriv);
-extern IndexingPrivate *asdf_prepare(IndexingMethod *indm,
- UnitCell *cell, struct detector *det,
- float *ltl);
+extern void *asdf_prepare(IndexingMethod *indm,
+ UnitCell *cell, struct detector *det,
+ float *ltl);
-extern void asdf_cleanup(IndexingPrivate *pp);
+extern void asdf_cleanup(void *pp);
#else /* HAVE_FFTW */
-int run_asdf(struct image *image, IndexingPrivate *ipriv)
+int run_asdf(struct image *image, void *ipriv)
{
ERROR("This copy of CrystFEL was compiled without FFTW support.\n");
return 0;
}
-IndexingPrivate *asdf_prepare(IndexingMethod *indm,
- UnitCell *cell, struct detector *det,
- float *ltl)
+void *asdf_prepare(IndexingMethod *indm,
+ UnitCell *cell, struct detector *det,
+ float *ltl)
{
ERROR("This copy of CrystFEL was compiled without FFTW support.\n");
ERROR("To use asdf indexing, recompile with FFTW.\n");
return NULL;
}
-void asdf_cleanup(IndexingPrivate *pp)
+void asdf_cleanup(void *pp)
{
}
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 8fa3b894..6610abc3 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -272,6 +272,12 @@ void cell_print(UnitCell *cell)
STATUS("c* = %10.3e %10.3e %10.3e m^-1 (modulus %10.3e m^-1)\n",
csx, csy, csz, modulus(csx, csy, csz));
+ STATUS("alpha* = %6.2f deg, beta* = %6.2f deg, "
+ "gamma* = %6.2f deg\n",
+ rad2deg(angle_between(bsx, bsy, bsz, csx, csy, csz)),
+ rad2deg(angle_between(asx, asy, asz, csx, csy, csz)),
+ rad2deg(angle_between(asx, asy, asz, bsx, bsy, bsz)));
+
STATUS("Cell representation is %s.\n", cell_rep(cell));
} else {
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index ec591e24..3de61073 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -460,10 +460,12 @@ int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c,
case CELL_REP_RECIP:
/* Convert reciprocal -> crystallographic.
* Start by converting reciprocal -> cartesian */
- cell_invert(cell->axs, cell->ays, cell->azs,
- cell->bxs, cell->bys, cell->bzs,
- cell->cxs, cell->cys, cell->czs,
- &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+ if ( cell_invert(cell->axs, cell->ays, cell->azs,
+ cell->bxs, cell->bys, cell->bzs,
+ cell->cxs, cell->cys, cell->czs,
+ &ax, &ay, &az,
+ &bx, &by, &bz,
+ &cx, &cy, &cz) ) return 1;
/* Now convert cartesian -> crystallographic */
*a = modulus(ax, ay, az);
diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
index 42ffbe40..e3b351ff 100644
--- a/libcrystfel/src/detector.c
+++ b/libcrystfel/src/detector.c
@@ -31,8 +31,6 @@
*
*/
-#define _ISOC99_SOURCE
-#define _GNU_SOURCE
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
@@ -530,37 +528,15 @@ int panel_number(struct detector *det, struct panel *p)
}
-void fill_in_values(struct detector *det, struct hdfile *f, struct event* ev)
+void adjust_centering_for_rail(struct panel *p)
{
- int i;
-
- for ( i=0; i<det->n_panels; i++ ) {
-
- double offs;
- struct panel *p = &det->panels[i];
-
- if ( p->clen_from != NULL ) {
-
- double val;
- int r;
+ double offs;
- 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 = val * 1.0e-3;
- }
-
- }
-
- /* Offset in +z direction from calibrated clen to actual */
- offs = p->clen - p->clen_for_centering;
- p->cnx += p->rail_x * offs;
- p->cny += p->rail_y * offs;
- p->clen = p->clen_for_centering + p->coffset + p->rail_z * offs;
-
- }
+ /* Offset in +z direction from calibrated clen to actual */
+ offs = p->clen - p->clen_for_centering;
+ p->cnx += p->rail_x * offs;
+ p->cny += p->rail_y * offs;
+ p->clen = p->clen_for_centering + p->coffset + p->rail_z * offs;
}
@@ -911,7 +887,6 @@ static int parse_field_for_panel(struct panel *panel, const char *key,
panel->adu_per_eV = atof(val);
} else if ( strcmp(key, "adu_per_photon") == 0 ) {
panel->adu_per_photon = atof(val);
- STATUS("got adu per photon: %s\n", val);
} else if ( strcmp(key, "rigid_group") == 0 ) {
add_to_rigid_group(find_or_add_rg(det, val), panel);
} else if ( strcmp(key, "clen") == 0 ) {
@@ -1079,12 +1054,15 @@ static void parse_toplevel(struct detector *det, struct beam_params *beam,
} else if ( strcmp(key, "photon_energy") == 0 ) {
if ( beam != NULL ) {
- if ( strncmp(val, "/", 1) == 0 ) {
+ double v;
+ char *end;
+ v = strtod(val, &end);
+ if ( (val[0] != '\0') && (end[0] == '\0') ) {
+ beam->photon_energy = v;
+ beam->photon_energy_from = NULL;
+ } else {
beam->photon_energy = 0.0;
beam->photon_energy_from = strdup(val);
- } else {
- beam->photon_energy = atof(val);
- beam->photon_energy_from = NULL;
}
}
@@ -1204,7 +1182,7 @@ struct detector *get_detector_geometry_2(const char *filename,
int i;
int rgi, rgci;
int reject = 0;
- int path_dim;
+ int path_dim, mask_path_dim;
int dim_dim;
int dim_reject = 0;
int dim_dim_reject = 0;
@@ -1389,6 +1367,7 @@ struct detector *get_detector_geometry_2(const char *filename,
}
+ mask_path_dim = -1;
for ( i=0; i<det->n_panels; i++ ) {
int panel_mask_dim = 0;
@@ -1398,8 +1377,7 @@ struct detector *get_detector_geometry_2(const char *filename,
next_instance = det->panels[i].mask;
- while(next_instance)
- {
+ while ( next_instance ) {
next_instance = strstr(next_instance, "%");
if ( next_instance != NULL ) {
next_instance += 1*sizeof(char);
@@ -1407,18 +1385,29 @@ struct detector *get_detector_geometry_2(const char *filename,
}
}
- if ( panel_mask_dim != path_dim ) {
- dim_reject = 1;
+ if ( mask_path_dim == -1 ) {
+ mask_path_dim = panel_mask_dim;
+ } else {
+ if ( panel_mask_dim != mask_path_dim ) {
+ dim_reject = 1;
+ }
}
+
}
}
- if ( dim_reject == 1) {
+ if ( dim_reject == 1 ) {
ERROR("All panels' data and mask entries must have the same "
"number of placeholders\n");
reject = 1;
}
+ if ( mask_path_dim > path_dim ) {
+ ERROR("Number of placeholders in mask cannot be larger than "
+ "for data\n");
+ reject = 1;
+ }
+
det->path_dim = path_dim;
dim_dim_reject = 0;
@@ -1993,6 +1982,11 @@ double largest_q(struct image *image)
struct rvec q;
double tt;
+ if ( image->det == NULL ) {
+ ERROR("No detector geometry. assuming detector is infinite!\n");
+ return INFINITY;
+ }
+
q = get_q_for_panel(image->det->furthest_out_panel,
image->det->furthest_out_fs,
image->det->furthest_out_ss,
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
index 04e3c1ec..a2be2b47 100644
--- a/libcrystfel/src/detector.h
+++ b/libcrystfel/src/detector.h
@@ -3,12 +3,12 @@
*
* Detector properties
*
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
* Authors:
- * 2009-2016 Thomas White <taw@physics.org>
+ * 2009-2017 Thomas White <taw@physics.org>
* 2011-2012 Richard Kirian <rkirian@asu.edu>
* 2014 Valerio Mariani
* 2011 Andrew Aquila
@@ -42,7 +42,6 @@ struct rg_collection;
struct detector;
struct panel;
struct badregion;
-struct detector;
struct beam_params;
struct hdfile;
struct event;
@@ -80,6 +79,47 @@ struct rg_collection
};
+/**
+ * panel:
+ * @name: Text name for the panel (fixed length array)
+ * @cnx: Location of corner, in pixels, x coordinate
+ * @cny: Location of corner, in pixels, y coordinate
+ * @coffset: The offset to be applied from @clen (which may come from elsewhere)
+ * @clen: The distance from the interaction point to the corner of the first pixel
+ * @clen_from: Location to get @clen from, e.g. from HDF5 file
+ * @mask: Location of mask data
+ * @mask_file: Filename for mask data
+ * @satmap: Location of per-pixel saturation map
+ * @satmap_file: Filename for saturation map
+ * @res: Resolution of panel in pixels per metre
+ * @badrow: Readout direction (for filtering out clusters of peaks)
+ * @no_index: Non-zero if panel is entirely "bad"
+ * @adu_per_photon: Number of detector intensity units per photon
+ * @adu_per_eV: Number of detector intensity units per eV of photon energy
+ * @max_adu: Saturation value
+ * @dim_structure: Dimension structure
+ * @fsx: Real-space x-direction of data fast-scan direction
+ * @fsy: Real-space y-direction of data fast-scan direction
+ * @fsz: Real-space z-direction of data fast-scan direction
+ * @ssx: Real-space x-direction of data slow-scan direction
+ * @ssy: Real-space y-direction of data slow-scan direction
+ * @ssz: Real-space z-direction of data slow-scan direction
+ * @rail_x: x direction of camera length "rail"
+ * @rail_y: y direction of camera length "rail"
+ * @rail_z: z direction of camera length "rail"
+ * @clen_for_centering: Value of clen (without coffset) at which beam is centered
+ * @xfs: Data fast-scan direction of real-space x-direction
+ * @yfs: Data fast-scan direction of real-space y-direction
+ * @xss: Data slow-scan direction of real-space x-direction
+ * @yss: Data slow-scan direction of real-space y-direction
+ * @orig_min_fs: Minimum fs coordinate of data in file
+ * @orig_max_fs: Maximum fs coordinate of data in file
+ * @orig_min_ss: Minimum ss coordinate of data in file (inclusive)
+ * @orig_max_ss: Maximum ss coordinate of data in file (inclusive)
+ * @data: Location of data in file
+ * @w: Width of panel
+ * @h: Height of panel
+ */
struct panel
{
char name[1024]; /* Name for this panel */
@@ -224,9 +264,8 @@ extern void get_pixel_extents(struct detector *det,
double *min_x, double *min_y,
double *max_x, double *max_y);
-extern void fill_in_values(struct detector *det, struct hdfile *f,
- struct event* ev);
extern void fill_in_adu(struct image *image);
+extern void adjust_centering_for_rail(struct panel *p);
extern int panel_is_in_rigid_group(const struct rigid_group *rg,
struct panel *p);
diff --git a/libcrystfel/src/dirax.c b/libcrystfel/src/dirax.c
index 19f35696..e9466a24 100644
--- a/libcrystfel/src/dirax.c
+++ b/libcrystfel/src/dirax.c
@@ -3,11 +3,11 @@
*
* Invoke the DirAx auto-indexing program
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2014 Thomas White <taw@physics.org>
+ * 2010-2014,2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -532,7 +532,7 @@ static void write_drx(struct image *image)
}
-int run_dirax(struct image *image, IndexingPrivate *ipriv)
+int run_dirax(struct image *image, void *ipriv)
{
unsigned int opts;
int status;
@@ -638,8 +638,8 @@ int run_dirax(struct image *image, IndexingPrivate *ipriv)
}
-IndexingPrivate *dirax_prepare(IndexingMethod *indm, UnitCell *cell,
- struct detector *det, float *ltl)
+void *dirax_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct detector *det, float *ltl)
{
struct dirax_private *dp;
int need_cell = 0;
@@ -670,7 +670,7 @@ IndexingPrivate *dirax_prepare(IndexingMethod *indm, UnitCell *cell,
}
-void dirax_cleanup(IndexingPrivate *pp)
+void dirax_cleanup(void *pp)
{
struct dirax_private *p;
p = (struct dirax_private *)pp;
diff --git a/libcrystfel/src/dirax.h b/libcrystfel/src/dirax.h
index 9776f3f0..96ba6dbc 100644
--- a/libcrystfel/src/dirax.h
+++ b/libcrystfel/src/dirax.h
@@ -3,11 +3,11 @@
*
* Invoke the DirAx auto-indexing program
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010,2012-2014 Thomas White <taw@physics.org>
+ * 2010,2012-2014,2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -39,13 +39,12 @@
extern "C" {
#endif
-extern int run_dirax(struct image *image, IndexingPrivate *ipriv);
+extern int run_dirax(struct image *image, void *ipriv);
-extern IndexingPrivate *dirax_prepare(IndexingMethod *indm,
- UnitCell *cell, struct detector *det,
- float *ltl);
+extern void *dirax_prepare(IndexingMethod *indm,
+ UnitCell *cell, struct detector *det, float *ltl);
-extern void dirax_cleanup(IndexingPrivate *pp);
+extern void dirax_cleanup(void *pp);
#ifdef __cplusplus
}
diff --git a/libcrystfel/src/events.c b/libcrystfel/src/events.c
index 25771a69..8e4eb861 100644
--- a/libcrystfel/src/events.c
+++ b/libcrystfel/src/events.c
@@ -3,10 +3,11 @@
*
* Event properties
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
+ * 2017 Thomas White
* 2014 Valerio Mariani
*
* This file is part of CrystFEL.
@@ -26,10 +27,8 @@
*
*/
-#define _ISOC99_SOURCE
-#define _GNU_SOURCE
-
#include "events.h"
+#include "utils.h"
#include <hdf5.h>
#include <string.h>
@@ -585,46 +584,38 @@ char *retrieve_full_path(struct event *ev, const char *data)
{
int ei ;
char *return_value;
+ char *pholder;
return_value = strdup(data);
+ pholder = strstr(return_value,"%");
+ ei = 0;
- for ( ei=0; ei<ev->path_length; ei++ ) {
+ while ( pholder != NULL ) {
char *tmp;
- tmp = event_path_placeholder_subst(ev->path_entries[ei],
- return_value);
- free(return_value);
- return_value = tmp;
-
- }
-
- return return_value;
-
-}
-
-
-char *partial_event_substitution(struct event *ev, const char *data)
-{
- int ei ;
- char *return_value;
- char *pholder;
-
- return_value = strdup(data);
- pholder = strstr(return_value,"%");
- ei = 0;
+ /* Check we have enough things to put in the placeholders */
+ if ( ei >= ev->path_length ) {
+ ERROR("Too many placeholders ('%%') in location.\n");
+ return NULL;
+ }
- while( pholder != NULL) {
+ /* Substitute one placeholder */
+ tmp = event_path_placeholder_subst(ev->path_entries[ei++],
+ return_value);
- char *tmp_subst_data;
+ if ( tmp == NULL ) {
+ ERROR("Couldn't substitute placeholder\n");
+ return NULL;
+ }
- tmp_subst_data = event_path_placeholder_subst(ev->path_entries[ei],
- return_value);
+ /* Next time, we will substitute the next part of the path into
+ * the partially substituted string */
free(return_value);
- return_value = strdup(tmp_subst_data);
- free(tmp_subst_data);
+ return_value = tmp;
+
pholder = strstr(return_value, "%");
- ei += 1;
+
}
return return_value;
diff --git a/libcrystfel/src/events.h b/libcrystfel/src/events.h
index 7f9c6731..4f717209 100644
--- a/libcrystfel/src/events.h
+++ b/libcrystfel/src/events.h
@@ -78,7 +78,6 @@ extern char *get_event_string(struct event *ev);
extern struct event *get_event_from_event_string(const char *ev_string);
extern char *event_path_placeholder_subst(const char *ev_name,
const char *data);
-extern char *partial_event_substitution(struct event *ev, const char *data);
extern char *retrieve_full_path(struct event *ev, const char *data);
diff --git a/libcrystfel/src/felix.c b/libcrystfel/src/felix.c
index 02d523c6..b531e3af 100644
--- a/libcrystfel/src/felix.c
+++ b/libcrystfel/src/felix.c
@@ -3,8 +3,8 @@
*
* Invoke Felix for multi-crystal autoindexing.
*
- * Copyright © 2015 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2015-2017 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
* 2015 Thomas White <taw@physics.org>
@@ -658,9 +658,8 @@ static void parse_options(const char *options, struct felix_private *gp)
free(option);
}
-IndexingPrivate *felix_prepare(IndexingMethod *indm, UnitCell *cell,
- struct detector *det, float *ltl,
- const char *options)
+void *felix_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct detector *det, float *ltl, const char *options)
{
struct felix_private *gp;
diff --git a/libcrystfel/src/felix.h b/libcrystfel/src/felix.h
index 40771933..c06e963a 100644
--- a/libcrystfel/src/felix.h
+++ b/libcrystfel/src/felix.h
@@ -3,12 +3,12 @@
*
* Invoke Felix for multi-crystal autoindexing
*
- * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2013 Thomas White <taw@physics.org>
- * 2013-2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ * 2010-2013,2017 Thomas White <taw@physics.org>
+ * 2013-2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
*
* This file is part of CrystFEL.
*
@@ -36,11 +36,9 @@
#include "cell.h"
-extern IndexingPrivate *felix_prepare(IndexingMethod *indm,
- UnitCell *cell,
- struct detector *det,
- float *ltl,
- const char *options);
+extern void *felix_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct detector *det, float *ltl,
+ const char *options);
extern void felix_cleanup(IndexingPrivate *pp);
diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c
index 6542b360..13339506 100644
--- a/libcrystfel/src/hdf5-file.c
+++ b/libcrystfel/src/hdf5-file.c
@@ -44,6 +44,19 @@
#include "utils.h"
+/**
+ * SECTION:hdf5-file
+ * @short_description: HDF5 file handling
+ * @title: HDF5 file handling
+ * @section_id:
+ * @see_also:
+ * @include: "hdf5-file.h"
+ * @Image:
+ *
+ * Routines for accessing HDF5 files.
+ **/
+
+
struct hdf5_write_location {
const char *location;
@@ -121,8 +134,7 @@ struct hdfile *hdfile_open(const char *filename)
}
-int hdfile_set_image(struct hdfile *f, const char *path,
- struct panel *p)
+int hdfile_set_image(struct hdfile *f, const char *path)
{
f->dh = H5Dopen2(f->fh, path, H5P_DEFAULT);
if ( f->dh < 0 ) {
@@ -320,9 +332,34 @@ static float *read_hdf5_data(struct hdfile *f, char *path, int line)
}
-/* 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)
+/**
+ * get_peaks_cxi_2:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ * @fpe: A %filename_plus_event structure specifying the event
+ * @half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates
+ *
+ * Get peaks from HDF5, in "CXI format" (as in "CXIDB"). The data should be in
+ * a set of arrays under @p. The number of peaks should be in a 1D array at
+ * @p/nPeaks. The fast-scan and slow-scan coordinates should be in 2D arrays at
+ * @p/peakXPosRaw and @p/peakYPosRaw respectively (sorry about the naming). The
+ * first dimension of these arrays should be the event number (as given by
+ * @fpe). The intensities are expected to be at @p/peakTotalIntensity in a
+ * similar 2D array.
+ *
+ * CrystFEL considers all peak locations to be distances from the corner of the
+ * detector panel, in pixel units, consistent with its description of detector
+ * geometry (see 'man crystfel_geometry'). The software which generates the
+ * CXI files, including Cheetah, may instead consider the peak locations to be
+ * pixel indices in the data array. In this case, the peak coordinates should
+ * have 0.5 added to them. This will be done if @half_pixel_shift is non-zero.
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p,
+ struct filename_plus_event *fpe, int half_pixel_shift)
{
char path_n[1024];
char path_x[1024];
@@ -338,6 +375,8 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
float *buf_y;
float *buf_i;
+ double peak_offset = half_pixel_shift ? 0.5 : 0.0;
+
if ( (fpe != NULL) && (fpe->ev != NULL)
&& (fpe->ev->dim_entries != NULL) )
{
@@ -375,8 +414,8 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
float fs, ss, val;
struct panel *p;
- fs = buf_x[pk];
- ss = buf_y[pk];
+ fs = buf_x[pk] + peak_offset;
+ ss = buf_y[pk] + peak_offset;
val = buf_i[pk];
p = find_orig_panel(image->det, fs, ss);
@@ -395,7 +434,53 @@ int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
}
-int get_peaks(struct image *image, struct hdfile *f, const char *p)
+/**
+ * get_peaks_cxi:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ * @fpe: A %filename_plus_event structure specifying the event
+ *
+ * This is a wrapper function to preserve API compatibility with older CrystFEL
+ * versions. Use get_peaks_cxi_2() instead.
+ *
+ * This function is equivalent to get_peaks_cxi_2(@image, @f, @p, @fpe, 1).
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
+ struct filename_plus_event *fpe)
+{
+ return get_peaks_cxi_2(image, f, p, fpe, 1);
+}
+
+
+/**
+ * get_peaks_2:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ * @half_pixel_shift: Non-zero if 0.5 should be added to all peak coordinates
+ *
+ * Get peaks from HDF5. The peak list should be located at @p in the HDF5 file,
+ * a 2D array where the first dimension equals the number of peaks and second
+ * dimension is three. The first two columns contain the fast scan and slow
+ * scan coordinates, respectively, of the peaks. The third column contains the
+ * intensities.
+ *
+ * CrystFEL considers all peak locations to be distances from the corner of the
+ * detector panel, in pixel units, consistent with its description of detector
+ * geometry (see 'man crystfel_geometry'). The software which generates the
+ * CXI files, including Cheetah, may instead consider the peak locations to be
+ * pixel indices in the data array. In this case, the peak coordinates should
+ * have 0.5 added to them. This will be done if @half_pixel_shift is non-zero.
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks_2(struct image *image, struct hdfile *f, const char *p,
+ int half_pixel_shift)
{
hid_t dh, sh;
hsize_t size[2];
@@ -405,6 +490,7 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p)
herr_t r;
int tw;
char *np;
+ double peak_offset = half_pixel_shift ? 0.5 : 0.0;
if ( image->event != NULL ) {
np = retrieve_full_path(image->event, p);
@@ -473,8 +559,8 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p)
float fs, ss, val;
struct panel *p;
- fs = buf[tw*i+0];
- ss = buf[tw*i+1];
+ fs = buf[tw*i+0] + peak_offset;
+ ss = buf[tw*i+1] + peak_offset;
val = buf[tw*i+2];
p = find_orig_panel(image->det, fs, ss);
@@ -499,6 +585,26 @@ int get_peaks(struct image *image, struct hdfile *f, const char *p)
}
+/**
+ * get_peaks:
+ * @image: An %image structure
+ * @f: An %hdfile structure
+ * @p: The HDF5 path to the peak data
+ *
+ * This is a wrapper function to preserve API compatibility with older CrystFEL
+ * versions. Use get_peaks_2() instead.
+ *
+ * This function is equivalent to get_peaks_2(@image, @f, @p, 1).
+ *
+ * Returns: non-zero on error, zero otherwise.
+ *
+ */
+int get_peaks(struct image *image, struct hdfile *f, const char *p)
+{
+ return get_peaks_2(image, f, p, 1);
+}
+
+
static void cleanup(hid_t fh)
{
int n_ids, i;
@@ -1127,6 +1233,9 @@ static int get_scalar_value(struct hdfile *f, const char *name, void *val,
return 1;
}
+ H5Tclose(type);
+ H5Dclose(dh);
+
return 0;
}
@@ -1154,7 +1263,7 @@ static int get_ev_based_value(struct hdfile *f, const char *name,
char *subst_name = NULL;
if ( ev->path_length != 0 ) {
- subst_name = partial_event_substitution(ev, name);
+ subst_name = retrieve_full_path(ev, name);
} else {
subst_name = strdup(name);
}
@@ -1284,8 +1393,10 @@ int hdfile_get_value(struct hdfile *f, const char *name,
}
-void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f,
- struct event *ev, struct image *image)
+static void hdfile_fill_in_beam_parameters(struct beam_params *beam,
+ struct hdfile *f,
+ struct event *ev,
+ struct image *image)
{
double eV;
@@ -1298,8 +1409,8 @@ void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f,
int r;
- r = hdfile_get_value(f, beam->photon_energy_from, ev, &eV,
- H5T_NATIVE_DOUBLE);
+ 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);
@@ -1311,6 +1422,36 @@ void fill_in_beam_parameters(struct beam_params *beam, struct hdfile *f,
}
+static void hdfile_fill_in_clen(struct detector *det, struct hdfile *f,
+ struct event *ev)
+{
+ int i;
+
+ for ( i=0; i<det->n_panels; i++ ) {
+
+ struct panel *p = &det->panels[i];
+
+ if ( p->clen_from != NULL ) {
+
+ double val;
+ int r;
+
+ 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 = val * 1.0e-3;
+ }
+
+ }
+
+ adjust_centering_for_rail(p);
+
+ }
+}
+
+
int hdf5_read(struct hdfile *f, struct image *image, const char *element,
int satcorr)
{
@@ -1326,7 +1467,7 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element,
if ( element == NULL ) {
fail = hdfile_set_first_image(f, "/");
} else {
- fail = hdfile_set_image(f, element, NULL);
+ fail = hdfile_set_image(f, element);
}
if ( fail ) {
@@ -1375,7 +1516,7 @@ int hdf5_read(struct hdfile *f, struct image *image, const char *element,
if ( image->beam != NULL ) {
- fill_in_beam_parameters(image->beam, f, NULL, image);
+ hdfile_fill_in_beam_parameters(image->beam, f, NULL, image);
if ( image->lambda > 1000 ) {
/* Error message covers a silly value in the beam file
@@ -1634,10 +1775,13 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
if ( !exists ) {
ERROR("Cannot find data for panel %s\n",
p->name);
+ free(image->dp);
+ free(image->bad);
+ free(image->sat);
return 1;
}
- fail = hdfile_set_image(f, panel_full_path, p);
+ fail = hdfile_set_image(f, panel_full_path);
free(panel_full_path);
@@ -1654,9 +1798,12 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
if ( !exists ) {
ERROR("Cannot find data for panel %s\n",
p->name);
+ free(image->dp);
+ free(image->bad);
+ free(image->sat);
return 1;
}
- fail = hdfile_set_image(f, p->data, p);
+ fail = hdfile_set_image(f, p->data);
}
@@ -1664,6 +1811,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
if ( fail ) {
ERROR("Couldn't select path for panel %s\n",
p->name);
+ free(image->dp);
+ free(image->bad);
+ free(image->sat);
return 1;
}
@@ -1673,6 +1823,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
f_count = malloc(hsd->num_dims*sizeof(hsize_t));
if ( (f_offset == NULL) || (f_count == NULL ) ) {
ERROR("Failed to allocate offset or count.\n");
+ free(image->dp);
+ free(image->bad);
+ free(image->sat);
return 1;
}
for ( hsi=0; hsi<hsd->num_dims; hsi++ ) {
@@ -1700,6 +1853,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
if ( check < 0 ) {
ERROR("Error selecting file dataspace for panel %s\n",
p->name);
+ free(image->dp);
+ free(image->bad);
+ free(image->sat);
return 1;
}
@@ -1713,6 +1869,9 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
ERROR("Failed to allocate panel %s\n", p->name);
free(f_offset);
free(f_count);
+ free(image->dp);
+ free(image->bad);
+ free(image->sat);
return 1;
}
for ( i=0; i<p->w*p->h; i++ ) image->sat[pi][i] = INFINITY;
@@ -1724,6 +1883,13 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
p->name);
free(f_offset);
free(f_count);
+ for ( i=0; i<=pi; i++ ) {
+ free(image->dp[i]);
+ free(image->sat[i]);
+ }
+ free(image->dp);
+ free(image->bad);
+ free(image->sat);
return 1;
}
@@ -1740,7 +1906,7 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
if ( load_satmap(f, ev, p, f_offset, f_count, hsd,
image->sat[pi]) )
{
- ERROR("Failed to laod sat map for panel %s\n",
+ ERROR("Failed to load sat map for panel %s\n",
p->name);
}
}
@@ -1753,13 +1919,13 @@ int hdf5_read2(struct hdfile *f, struct image *image, struct event *ev,
H5Dclose(f->dh);
f->data_open = 0;
- fill_in_values(image->det, f, ev);
+ hdfile_fill_in_clen(image->det, f, ev);
if ( satcorr ) debodge_saturation(f, image);
if ( image->beam != NULL ) {
- fill_in_beam_parameters(image->beam, f, ev, image);
+ hdfile_fill_in_beam_parameters(image->beam, f, ev, image);
if ( (image->lambda > 1.0) || (image->lambda < 1e-20) ) {
@@ -1966,7 +2132,7 @@ char *hdfile_get_string_value(struct hdfile *f, const char *name,
char *tmp = NULL, *subst_name = NULL;
if (ev != NULL && ev->path_length != 0 ) {
- subst_name = partial_event_substitution(ev, name);
+ subst_name = retrieve_full_path(ev, name);
} else {
subst_name = strdup(name);
}
@@ -2166,7 +2332,7 @@ int hdfile_set_first_image(struct hdfile *f, const char *group)
for ( i=0; i<n; i++ ) {
if ( is_image[i] ) {
- hdfile_set_image(f, names[i], NULL);
+ hdfile_set_image(f, names[i]);
for ( j=0; j<n; j++ ) free(names[j]);
free(is_image);
free(is_group);
@@ -2484,6 +2650,9 @@ static int check_dims(struct hdfile *hdfile, struct panel *p, struct event *ev,
return 1;
}
+ H5Sclose(sh);
+ H5Dclose(dh);
+
return 0;
}
diff --git a/libcrystfel/src/hdf5-file.h b/libcrystfel/src/hdf5-file.h
index 8c89eb93..ab53dd2e 100644
--- a/libcrystfel/src/hdf5-file.h
+++ b/libcrystfel/src/hdf5-file.h
@@ -3,11 +3,11 @@
*
* Read/write HDF5 data files
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2009-2015 Thomas White <taw@physics.org>
+ * 2009-2017 Thomas White <taw@physics.org>
* 2014 Valerio Mariani
*
@@ -67,8 +67,8 @@ extern int hdf5_read2(struct hdfile *f, struct image *image,
extern int check_path_existence(hid_t fh, const char *path);
extern struct hdfile *hdfile_open(const char *filename);
-int hdfile_set_image(struct hdfile *f, const char *path,
- struct panel *p);
+int hdfile_set_image(struct hdfile *f, const char *path);
+
extern int16_t *hdfile_get_image_binned(struct hdfile *hdfile,
int binning, int16_t *maxp);
extern char **hdfile_read_group(struct hdfile *f, int *n, const char *parent,
@@ -78,9 +78,16 @@ extern void hdfile_close(struct hdfile *f);
extern int get_peaks(struct image *image, struct hdfile *f, const char *p);
+extern int get_peaks_2(struct image *image, struct hdfile *f, const char *p,
+ int half_pixel_shift);
+
extern int get_peaks_cxi(struct image *image, struct hdfile *f, const char *p,
struct filename_plus_event *fpe);
+extern int get_peaks_cxi_2(struct image *image, struct hdfile *f, const char *p,
+ struct filename_plus_event *fpe,
+ int half_pixel_shift);
+
extern struct copy_hdf5_field *new_copy_hdf5_field_list(void);
extern void free_copy_hdf5_field_list(struct copy_hdf5_field *f);
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
index 7bb4cede..aad5017c 100644
--- a/libcrystfel/src/image.c
+++ b/libcrystfel/src/image.c
@@ -3,12 +3,12 @@
*
* Handle images and image features
*
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
* 2014 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
- * 2011-2016 Thomas White <taw@physics.org>
+ * 2011-2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -27,14 +27,28 @@
*
*/
-#define _ISOC99_SOURCE
+#include <config.h>
+
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <stdio.h>
+#include <hdf5.h>
+
+#ifdef HAVE_CBFLIB
+#ifdef HAVE_CBFLIB_CBF_H
+#include <cbflib/cbf.h>
+#endif
+#ifdef HAVE_CBF_CBF_H
+#include <cbf/cbf.h>
+#endif
+#endif
#include "image.h"
#include "utils.h"
+#include "events.h"
+#include "hdf5-file.h"
+#include "detector.h"
/**
* SECTION:image
@@ -51,6 +65,14 @@
*/
+struct imagefile
+{
+ enum imagefile_type type;
+ char *filename;
+ struct hdfile *hdfile;
+};
+
+
struct _imagefeaturelist
{
struct imagefeature *features;
@@ -300,3 +322,673 @@ void free_all_crystals(struct image *image)
free(image->crystals);
image->n_crystals = 0;
}
+
+
+/**************************** Image field lists *******************************/
+
+struct imagefile_field_list
+{
+ char **fields;
+ int n_fields;
+ int max_fields;
+};
+
+
+struct imagefile_field_list *new_imagefile_field_list()
+{
+ struct imagefile_field_list *n;
+
+ n = calloc(1, sizeof(struct imagefile_field_list));
+ if ( n == NULL ) return NULL;
+
+ n->max_fields = 32;
+ n->fields = malloc(n->max_fields*sizeof(char *));
+ if ( n->fields == NULL ) {
+ free(n);
+ return NULL;
+ }
+
+ return n;
+}
+
+
+void free_imagefile_field_list(struct imagefile_field_list *n)
+{
+ int i;
+ for ( i=0; i<n->n_fields; i++ ) {
+ free(n->fields[i]);
+ }
+ free(n->fields);
+ free(n);
+}
+
+
+void add_imagefile_field(struct imagefile_field_list *copyme, const char *name)
+{
+ int i;
+
+ /* Already on the list? Don't re-add if so. */
+ for ( i=0; i<copyme->n_fields; i++ ) {
+ if ( strcmp(copyme->fields[i], name) == 0 ) return;
+ }
+
+ /* Need more space? */
+ if ( copyme->n_fields == copyme->max_fields ) {
+
+ char **nfields;
+ int nmax = copyme->max_fields + 32;
+
+ nfields = realloc(copyme->fields, nmax*sizeof(char *));
+ if ( nfields == NULL ) {
+ ERROR("Failed to allocate space for new HDF5 field.\n");
+ return;
+ }
+
+ copyme->max_fields = nmax;
+ copyme->fields = nfields;
+
+ }
+
+ copyme->fields[copyme->n_fields] = strdup(name);
+ if ( copyme->fields[copyme->n_fields] == NULL ) {
+ ERROR("Failed to add field for copying '%s'\n", name);
+ return;
+ }
+
+ copyme->n_fields++;
+}
+
+
+/******************************* CBF files ************************************/
+
+#ifdef HAVE_CBFLIB
+
+static char *cbf_strerr(int e)
+{
+ char *err;
+
+ err = malloc(1024);
+ if ( err == NULL ) return NULL;
+
+ err[0] = '\0';
+
+ /* NB Sum of lengths of all strings must be less than 1024 */
+ if ( e & CBF_FORMAT ) strcat(err, "Invalid format");
+ if ( e & CBF_ALLOC ) strcat(err, "Memory allocation failed");
+ if ( e & CBF_ARGUMENT ) strcat(err, "Invalid argument");
+ if ( e & CBF_ASCII ) strcat(err, "Value is ASCII");
+ if ( e & CBF_BINARY ) strcat(err, "Value is binary");
+ if ( e & CBF_BITCOUNT ) strcat(err, "Wrong number of bits");
+ if ( e & CBF_ENDOFDATA ) strcat(err, "End of data");
+ if ( e & CBF_FILECLOSE ) strcat(err, "File close error");
+ if ( e & CBF_FILEOPEN ) strcat(err, "File open error");
+ if ( e & CBF_FILEREAD ) strcat(err, "File read error");
+ if ( e & CBF_FILETELL ) strcat(err, "File tell error");
+ if ( e & CBF_FILEWRITE ) strcat(err, "File write error");
+ if ( e & CBF_IDENTICAL ) strcat(err, "Name already exists");
+ if ( e & CBF_NOTFOUND ) strcat(err, "Not found");
+ if ( e & CBF_OVERFLOW ) strcat(err, "Overflow");
+ if ( e & CBF_UNDEFINED ) strcat(err, "Number undefined");
+ if ( e & CBF_NOTIMPLEMENTED ) strcat(err, "Not implemented");
+
+ return err;
+}
+
+
+static int unpack_panels(struct image *image, signed int *data, int data_width,
+ int data_height)
+{
+ int pi;
+
+ /* FIXME: Load these masks from an HDF5 file, if filenames are
+ * given in the geometry file */
+ uint16_t *flags = NULL;
+ float *sat = NULL;
+
+ image->dp = malloc(image->det->n_panels * sizeof(float *));
+ image->bad = malloc(image->det->n_panels * sizeof(int *));
+ image->sat = malloc(image->det->n_panels * sizeof(float *));
+ if ( (image->dp == NULL) || (image->bad == NULL)
+ || (image->sat == NULL) )
+ {
+ ERROR("Failed to allocate panels.\n");
+ return 1;
+ }
+
+ for ( pi=0; pi<image->det->n_panels; pi++ ) {
+
+ struct panel *p;
+ int fs, ss;
+
+ p = &image->det->panels[pi];
+ image->dp[pi] = malloc(p->w*p->h*sizeof(float));
+ image->bad[pi] = calloc(p->w*p->h, sizeof(int));
+ image->sat[pi] = malloc(p->w*p->h*sizeof(float));
+ if ( (image->dp[pi] == NULL) || (image->bad[pi] == NULL)
+ || (image->sat[pi] == NULL) )
+ {
+ ERROR("Failed to allocate panel\n");
+ return 1;
+ }
+
+ if ( p->mask != NULL ) {
+ ERROR("WARNING: Bad pixel masks do not currently work "
+ "with CBF files\n");
+ ERROR(" (bad pixel regions specified in the geometry "
+ "file will be used, however)\n");
+ }
+
+ if ( p->satmap != NULL ) {
+ ERROR("WARNING: Saturation maps do not currently work "
+ "with CBF files\n");
+ }
+
+ if ( (p->orig_min_fs + p->w > data_width)
+ || (p->orig_min_ss + p->h > data_height) )
+ {
+ ERROR("Panel %s is outside range of data in CBF file\n",
+ p->name);
+ return 1;
+ }
+
+ for ( ss=0; ss<p->h; ss++ ) {
+ for ( fs=0; fs<p->w; fs++ ) {
+
+ int idx;
+ int cfs, css;
+ int bad = 0;
+
+ cfs = fs+p->orig_min_fs;
+ css = ss+p->orig_min_ss;
+ idx = cfs + css*data_width;
+
+ image->dp[pi][fs+p->w*ss] = data[idx];
+
+ if ( sat != NULL ) {
+ image->sat[pi][fs+p->w*ss] = sat[idx];
+ } else {
+ image->sat[pi][fs+p->w*ss] = INFINITY;
+ }
+
+ if ( p->no_index ) bad = 1;
+
+ if ( in_bad_region(image->det, p, cfs, css) ) {
+ bad = 1;
+ }
+
+ if ( flags != NULL ) {
+
+ int f;
+
+ f = flags[idx];
+
+ /* Bad if it's missing any of the "good" bits */
+ if ( (f & image->det->mask_good)
+ != image->det->mask_good ) bad = 1;
+
+ /* Bad if it has any of the "bad" bits. */
+ if ( f & image->det->mask_bad ) bad = 1;
+
+ }
+ image->bad[pi][fs+p->w*ss] = bad;
+
+ }
+ }
+
+ }
+
+ return 0;
+}
+
+
+static void cbf_fill_in_beam_parameters(struct beam_params *beam,
+ struct imagefile *f,
+ struct image *image)
+{
+ double eV;
+
+ if ( beam->photon_energy_from == NULL ) {
+
+ /* Explicit value given */
+ eV = beam->photon_energy;
+
+ } else {
+
+ ERROR("Can't get photon energy from CBF yet.\n");
+ eV = 0.0;
+
+ }
+
+ image->lambda = ph_en_to_lambda(eV_to_J(eV))*beam->photon_energy_scale;
+}
+
+
+static void cbf_fill_in_clen(struct detector *det, struct imagefile *f)
+{
+ int i;
+
+ for ( i=0; i<det->n_panels; i++ ) {
+
+ struct panel *p = &det->panels[i];
+
+ if ( p->clen_from != NULL ) {
+
+ ERROR("Can't get clen from CBF yet.\n");
+
+ }
+
+ adjust_centering_for_rail(p);
+
+ }
+}
+
+
+static int read_cbf(struct imagefile *f, struct image *image)
+{
+ cbf_handle cbfh;
+ FILE *fh;
+ int r;
+ unsigned int compression;
+ int binary_id, minelement, maxelement, elsigned, elunsigned;
+ size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding;
+ const char *byteorder;
+ signed int *data;
+
+ if ( image->det == NULL ) {
+ ERROR("read_cbf() needs a geometry\n");
+ return 1;
+ }
+
+ if ( cbf_make_handle(&cbfh) ) {
+ ERROR("Failed to allocate CBF handle\n");
+ return 1;
+ }
+
+ fh = fopen(f->filename, "rb");
+ if ( fh == NULL ) {
+ ERROR("Failed to open '%s'\n", f->filename);
+ return 1;
+ }
+ /* CBFlib calls fclose(fh) when it's ready */
+
+ if ( cbf_read_widefile(cbfh, fh, 0) ) {
+ ERROR("Failed to read CBF file '%s'\n", f->filename);
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ /* Select row 0 in data column inside array_data */
+ cbf_find_category(cbfh, "array_data");
+ cbf_find_column(cbfh, "data");
+ cbf_select_row(cbfh, 0);
+
+ /* Get parameters for array read */
+ r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id,
+ &elsize, &elsigned, &elunsigned,
+ &elements,
+ &minelement, &maxelement,
+ &byteorder,
+ &dimfast, &dimmid, &dimslow,
+ &padding);
+ if ( r ) {
+ char *err = cbf_strerr(r);
+ ERROR("Failed to read CBF array parameters: %s\n", err);
+ free(err);
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ if ( dimslow != 0 ) {
+ ERROR("CBF data array is 3D - don't know what to do with it\n");
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ if ( dimfast*dimmid*elsize > 10e9 ) {
+ ERROR("CBF data is far too big (%i x %i x %i bytes).\n",
+ (int)dimfast, (int)dimmid, (int)elsize);
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ if ( elsize != 4 ) {
+ STATUS("Don't know what to do with element size %i\n",
+ (int)elsize);
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ if ( strcmp(byteorder, "little_endian") != 0 ) {
+ STATUS("Don't know what to do with non-little-endian datan\n");
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ data = malloc(elsize*dimfast*dimmid);
+ if ( data == NULL ) {
+ ERROR("Failed to allocate memory for CBF data\n");
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1,
+ elsize*dimfast*dimmid, &elread);
+ if ( r ) {
+ char *err = cbf_strerr(r);
+ ERROR("Failed to read CBF array: %s\n", err);
+ free(err);
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ unpack_panels(image, data, dimfast, dimmid);
+ free(data);
+
+ if ( image->beam != NULL ) {
+ cbf_fill_in_beam_parameters(image->beam, f, image);
+ if ( image->lambda > 1000 ) {
+ ERROR("WARNING: Missing or nonsensical wavelength "
+ "(%e m) for %s.\n",
+ image->lambda, image->filename);
+ }
+ }
+ cbf_fill_in_clen(image->det, f);
+ fill_in_adu(image);
+
+ cbf_free_handle(cbfh);
+ return 0;
+}
+
+
+static float *convert_float(signed int *data, int w, int h)
+{
+ float *df;
+ long int i;
+
+ df = malloc(sizeof(float)*w*h);
+ if ( df == NULL ) return NULL;
+
+ for ( i=0; i<w*h; i++ ) {
+ df[i] = data[i];
+ }
+
+ return df;
+}
+
+
+static int read_cbf_simple(struct imagefile *f, struct image *image)
+{
+ cbf_handle cbfh;
+ FILE *fh;
+ int r;
+ unsigned int compression;
+ int binary_id, minelement, maxelement, elsigned, elunsigned;
+ size_t elsize, elements, elread, dimfast, dimmid, dimslow, padding;
+ const char *byteorder;
+ signed int *data;
+
+ if ( cbf_make_handle(&cbfh) ) {
+ ERROR("Failed to allocate CBF handle\n");
+ return 1;
+ }
+
+ fh = fopen(f->filename, "rb");
+ if ( fh == NULL ) {
+ ERROR("Failed to open '%s'\n", f->filename);
+ return 1;
+ }
+ /* CBFlib calls fclose(fh) when it's ready */
+
+ if ( cbf_read_widefile(cbfh, fh, 0) ) {
+ ERROR("Failed to read CBF file '%s'\n", f->filename);
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ /* Select row 0 in data column inside array_data */
+ cbf_find_category(cbfh, "array_data");
+ cbf_find_column(cbfh, "data");
+ cbf_select_row(cbfh, 0);
+
+ /* Get parameters for array read */
+ r = cbf_get_integerarrayparameters_wdims(cbfh, &compression, &binary_id,
+ &elsize, &elsigned, &elunsigned,
+ &elements,
+ &minelement, &maxelement,
+ &byteorder,
+ &dimfast, &dimmid, &dimslow,
+ &padding);
+ if ( r ) {
+ char *err = cbf_strerr(r);
+ ERROR("Failed to read CBF array parameters: %s\n", err);
+ free(err);
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ if ( dimslow != 0 ) {
+ ERROR("CBF data array is 3D - don't know what to do with it\n");
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ if ( dimfast*dimmid*elsize > 10e9 ) {
+ ERROR("CBF data is far too big (%i x %i x %i bytes).\n",
+ (int)dimfast, (int)dimmid, (int)elsize);
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ if ( elsize != 4 ) {
+ STATUS("Don't know what to do with element size %i\n",
+ (int)elsize);
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ if ( strcmp(byteorder, "little_endian") != 0 ) {
+ STATUS("Don't know what to do with non-little-endian datan\n");
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ data = malloc(elsize*dimfast*dimmid);
+ if ( data == NULL ) {
+ ERROR("Failed to allocate memory for CBF data\n");
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ r = cbf_get_integerarray(cbfh, &binary_id, data, elsize, 1,
+ elsize*dimfast*dimmid, &elread);
+ if ( r ) {
+ char *err = cbf_strerr(r);
+ ERROR("Failed to read CBF array: %s\n", err);
+ free(err);
+ cbf_free_handle(cbfh);
+ return 1;
+ }
+
+ image->det = simple_geometry(image, dimfast, dimmid);
+ image->dp = malloc(sizeof(float *));
+ if ( image->dp == NULL ) {
+ ERROR("Failed to allocate dp array\n");
+ return 1;
+ }
+ image->dp[0] = convert_float(data, dimfast, dimmid);
+ if ( image->dp[0] == NULL ) {
+ ERROR("Failed to allocate dp array\n");
+ return 1;
+ }
+
+ if ( image->beam != NULL ) {
+ cbf_fill_in_beam_parameters(image->beam, f, image);
+ if ( image->lambda > 1000 ) {
+ ERROR("WARNING: Missing or nonsensical wavelength "
+ "(%e m) for %s.\n",
+ image->lambda, image->filename);
+ }
+ }
+ cbf_fill_in_clen(image->det, f);
+ fill_in_adu(image);
+
+ cbf_free_handle(cbfh);
+ return 0;
+}
+
+#else /* HAVE_CBFLIB */
+
+static int read_cbf_simple(struct imagefile *f, struct image *image)
+{
+ ERROR("This version of CrystFEL was compiled without CBF support.\n");
+ return 1;
+}
+
+
+static int read_cbf(struct imagefile *f, struct image *image)
+{
+ ERROR("This version of CrystFEL was compiled without CBF support.\n");
+ return 1;
+}
+
+
+#endif /* HAVE_CBFLIB */
+
+
+/****************************** Image files ***********************************/
+
+
+static signed int is_cbf_file(const char *filename)
+{
+ FILE *fh;
+ char line[1024];
+
+ fh = fopen(filename, "r");
+ if ( fh == NULL ) return -1;
+
+ if ( fgets(line, 1024, fh) == NULL ) return -1;
+ fclose(fh);
+
+ if ( strstr(line, "CBF") == NULL ) {
+ return 0;
+ }
+
+ return 1;
+}
+
+
+struct imagefile *imagefile_open(const char *filename)
+{
+ struct imagefile *f;
+
+ f = malloc(sizeof(struct imagefile));
+ if ( f == NULL ) return NULL;
+
+ if ( H5Fis_hdf5(filename) > 0 ) {
+
+ /* This is an HDF5, pass through to HDF5 layer */
+ f->type = IMAGEFILE_HDF5;
+ f->hdfile = hdfile_open(filename);
+
+ if ( f->hdfile == NULL ) {
+ free(f);
+ return NULL;
+ }
+
+ } else if ( is_cbf_file(filename) > 0 ) {
+
+ f->type = IMAGEFILE_CBF;
+
+ }
+
+ f->filename = strdup(filename);
+ return f;
+}
+
+
+int imagefile_read(struct imagefile *f, struct image *image,
+ struct event *event)
+{
+ if ( f->type == IMAGEFILE_HDF5 ) {
+ return hdf5_read2(f->hdfile, image, event, 0);
+ } else if ( f->type == IMAGEFILE_CBF ) {
+ return read_cbf(f, image);
+ } else {
+ ERROR("Unknown file type %i\n", f->type);
+ return 1;
+ }
+}
+
+
+/* Read a simple file, no multi-event, no prior geometry etc, and
+ * generate a geometry for it */
+int imagefile_read_simple(struct imagefile *f, struct image *image)
+{
+ if ( f->type == IMAGEFILE_HDF5 ) {
+ return hdf5_read(f->hdfile, image, NULL, 0);
+ } else {
+ return read_cbf_simple(f, image);
+ }
+}
+
+
+enum imagefile_type imagefile_get_type(struct imagefile *f)
+{
+ assert(f != NULL);
+ return f->type;
+}
+
+
+struct hdfile *imagefile_get_hdfile(struct imagefile *f)
+{
+ if ( f->type != IMAGEFILE_HDF5 ) {
+ ERROR("Not an HDF5 file!\n");
+ return NULL;
+ }
+
+ return f->hdfile;
+}
+
+
+void imagefile_copy_fields(struct imagefile *f,
+ const struct imagefile_field_list *copyme,
+ FILE *fh, struct event *ev)
+{
+ int i;
+
+ if ( copyme == NULL ) return;
+
+ for ( i=0; i<copyme->n_fields; i++ ) {
+
+ char *val;
+ char *field;
+
+ field = copyme->fields[i];
+
+ if ( f->type == IMAGEFILE_HDF5 ) {
+ val = hdfile_get_string_value(f->hdfile, field, ev);
+ if ( field[0] == '/' ) {
+ fprintf(fh, "hdf5%s = %s\n", field, val);
+ } else {
+ fprintf(fh, "hdf5/%s = %s\n", field, val);
+ }
+ free(val);
+
+ } else {
+ STATUS("Mock CBF variable\n");
+ fprintf(fh, "cbf/%s = %s\n", field, "(FIXME)");
+ }
+
+
+ }
+}
+
+
+void imagefile_close(struct imagefile *f)
+{
+ if ( f->type == IMAGEFILE_HDF5 ) {
+ hdfile_close(f->hdfile);
+ }
+ free(f->filename);
+ free(f);
+}
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 9fd9b495..9719bb59 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -3,11 +3,11 @@
*
* Handle images and image features
*
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2009-2016 Thomas White <taw@physics.org>
+ * 2009-2017 Thomas White <taw@physics.org>
* 2014 Valerio Mariani
*
*
@@ -44,6 +44,8 @@ struct detector;
struct imagefeature;
struct sample;
struct image;
+struct imagefile;
+struct imagefile_field_list;
#include "utils.h"
#include "cell.h"
@@ -51,6 +53,7 @@ struct image;
#include "reflist.h"
#include "crystal.h"
#include "index.h"
+#include "events.h"
/**
* SpectrumType:
@@ -67,7 +70,26 @@ typedef enum {
SPECTRUM_TWOCOLOUR
} SpectrumType;
-/* Structure describing a feature in an image */
+
+/**
+ * imagefeature:
+ * @parent: Image this feature belongs to
+ * @fs: Fast scan coordinate
+ * @ss: Slow scan coordinate
+ * @p: Pointer to panel
+ * @intensity: Intensity of peak
+ * @rx: Reciprocal x coordinate in m^-1
+ * @ry: Reciprocal y coordinate in m^-1
+ * @rz: Reciprocal z coordinate in m^-1
+ * @name: Text name for feature
+ *
+ * Represents a peak in an image.
+ *
+ * Note carefully that the @fs and @ss coordinates are the distances, measured
+ * in pixels, from the corner of the panel. They are NOT pixel indices.
+ * If the peak is in the middle of the first pixel, its coordinates would be
+ * 0.5,0.5.
+ */
struct imagefeature {
struct image *parent;
@@ -81,12 +103,21 @@ struct imagefeature {
double ry;
double rz;
- /* Internal use only */
+ const char *name;
+
+ /*< private >*/
int valid;
+};
- const char *name;
+
+/* An enum representing the image file formats we can handle */
+enum imagefile_type
+{
+ IMAGEFILE_HDF5,
+ IMAGEFILE_CBF
};
+
/* An opaque type representing a list of image features */
typedef struct _imagefeaturelist ImageFeatureList;
@@ -98,44 +129,46 @@ struct sample
};
+/**
+ * beam_params:
+ * @photon_energy: eV per photon
+ * @photon_energy_from: HDF5 dataset name
+ * @photon_energy_scale: Scale factor for photon energy, if it comes from HDF5
+ */
struct beam_params
{
- double photon_energy; /* eV per photon */
- char *photon_energy_from; /* HDF5 dataset name */
- double photon_energy_scale; /* Scale factor for photon energy, if the
- * energy is to be from the HDF5 file */
+ double photon_energy;
+ char *photon_energy_from;
+ double photon_energy_scale;
};
/**
* image:
- *
- * <programlisting>
- * struct image
- * {
- * Crystal **crystals;
- * int n_crystals;
- * IndexingMethod indexed_by;
- *
- * struct detector *det;
- * struct beam_params *beam;
- * char *filename;
- * const struct copy_hdf5_field *copyme;
- *
- * int id;
- *
- * double lambda;
- * double div;
- * double bw;
- *
- * int width;
- * int height;
- *
- * long long int num_peaks;
- * long long int num_saturated_peaks;
- * ImageFeatureList *features;
- * };
- * </programlisting>
+ * @crystals: Array of crystals in the image
+ * @n_crystals: The number of crystals in the image
+ * @indexed_by: Indexing method which indexed this pattern
+ * @det: Detector structure
+ * @beam: Beam parameters structure
+ * @filename: Filename for the image file
+ * @copyme: Fields to copy from the image file to the stream
+ * @id: ID number of the thread handling this image
+ * @serial: Serial number for this image
+ * @lambda: Wavelength
+ * @div: Divergence
+ * @bw: Bandwidth
+ * @num_peaks: The number of peaks
+ * @num_saturated_peaks: The number of saturated peaks
+ * @features: The peaks found in the image
+ * @dp: The image data, by panel
+ * @bad: The bad pixel mask, array by panel
+ * @sat: The per-pixel saturation mask, array by panel
+ * @event: Event ID for the image
+ * @stuff_from_stream: Items read back from the stream
+ * @avg_clen: Mean of camera length values for all panels
+ * @spectrum: Spectrum information
+ * @nsamples: Number of spectrum samples
+ * @spectrum_size: SIze of spectrum array
*
* The field <structfield>data</structfield> contains the raw image data, if it
* is currently available. The data might be available throughout the
@@ -152,8 +185,8 @@ struct beam_params
* returned by the low-level indexing system. <structfield>n_crystals</structfield>
* is the number of crystals which were found in the image.
*
- * <structfield>copyme</structfield> represents a list of HDF5 fields to copy
- * to the output stream.
+ * <structfield>copyme</structfield> represents a list of fields in the image
+ * file (e.g. HDF5 fields or CBF headers) to copy to the output stream.
**/
struct image;
@@ -171,7 +204,7 @@ struct image {
struct beam_params *beam; /* The nominal beam parameters */
char *filename;
struct event *event;
- const struct copy_hdf5_field *copyme;
+ const struct imagefile_field_list *copyme;
struct stuff_from_stream *stuff_from_stream;
double avg_clen; /* Average camera length extracted
@@ -192,8 +225,8 @@ struct image {
double bw; /* Bandwidth as a fraction */
/* Detected peaks */
- long long int num_peaks;
- long long int num_saturated_peaks;
+ long long num_peaks;
+ long long num_saturated_peaks;
ImageFeatureList *features;
};
@@ -233,6 +266,25 @@ extern void image_add_crystal(struct image *image, Crystal *cryst);
extern void remove_flagged_crystals(struct image *image);
extern void free_all_crystals(struct image *image);
+/* Image files */
+extern struct imagefile *imagefile_open(const char *filename);
+extern int imagefile_read(struct imagefile *f, struct image *image,
+ struct event *event);
+extern int imagefile_read_simple(struct imagefile *f, struct image *image);
+extern struct hdfile *imagefile_get_hdfile(struct imagefile *f);
+extern enum imagefile_type imagefile_get_type(struct imagefile *f);
+extern void imagefile_copy_fields(struct imagefile *f,
+ const struct imagefile_field_list *copyme,
+ FILE *fh, struct event *ev);
+extern void imagefile_close(struct imagefile *f);
+
+/* Field lists */
+extern struct imagefile_field_list *new_imagefile_field_list(void);
+extern void free_imagefile_field_list(struct imagefile_field_list *f);
+
+extern void add_imagefile_field(struct imagefile_field_list *copyme,
+ const char *name);
+
#ifdef __cplusplus
}
#endif
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index e7fd879d..bdabb062 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -59,6 +59,17 @@
#include "taketwo.h"
+struct _indexingprivate
+{
+ UnitCell *target_cell;
+ float tolerance[4];
+
+ int n_methods;
+ IndexingMethod *methods;
+ void **engine_private;
+};
+
+
static int debug_index(struct image *image)
{
Crystal *cr = crystal_new();
@@ -72,78 +83,119 @@ static int debug_index(struct image *image)
}
-IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
- struct detector *det, float *ltl,
- const char *options)
+static void *prepare_method(IndexingMethod *m, UnitCell *cell,
+ struct detector *det, float *ltl,
+ const char *options)
{
- int n;
- int nm = 0;
- IndexingPrivate **iprivs;
+ char *str;
+ IndexingMethod in = *m;
+ void *priv = NULL;
- while ( indm[nm] != INDEXING_NONE ) nm++;
- iprivs = malloc((nm+1) * sizeof(IndexingPrivate *));
+ switch ( *m & INDEXING_METHOD_MASK ) {
- for ( n=0; n<nm; n++ ) {
+ case INDEXING_NONE :
+ priv = "none";
+ break;
- int i;
- IndexingMethod in;
- char *str;
+ case INDEXING_DIRAX :
+ priv = dirax_prepare(m, cell, det, ltl);
+ break;
- in = indm[n];
+ case INDEXING_ASDF :
+ priv = asdf_prepare(m, cell, det, ltl);
+ break;
- switch ( indm[n] & INDEXING_METHOD_MASK ) {
+ case INDEXING_MOSFLM :
+ priv = mosflm_prepare(m, cell, det, ltl);
+ break;
- case INDEXING_DIRAX :
- iprivs[n] = dirax_prepare(&indm[n], cell, det, ltl);
- break;
+ case INDEXING_XDS :
+ priv = xds_prepare(m, cell, det, ltl);
+ break;
- case INDEXING_ASDF :
- iprivs[n] = asdf_prepare(&indm[n], cell, det, ltl);
- break;
+ case INDEXING_DEBUG :
+ priv = (IndexingPrivate *)strdup("Hello!");
+ break;
- case INDEXING_MOSFLM :
- iprivs[n] = mosflm_prepare(&indm[n], cell, det, ltl);
- break;
+ case INDEXING_FELIX :
+ priv = felix_prepare(m, cell, det, ltl, options);
+ break;
- case INDEXING_XDS :
- iprivs[n] = xds_prepare(&indm[n], cell, det, ltl);
- break;
+ case INDEXING_TAKETWO :
+ priv = taketwo_prepare(m, 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", *m);
+ break;
- case INDEXING_FELIX :
- iprivs[n] = felix_prepare(&indm[n], cell, det, ltl,
- options);
- break;
+ }
- case INDEXING_TAKETWO :
- iprivs[n] = taketwo_prepare(&indm[n], cell, det, ltl);
- break;
+ str = indexer_str(*m);
- default :
- ERROR("Don't know how to prepare indexing method %i\n",
- indm[n]);
- break;
+ if ( priv == NULL ) {
+ ERROR("Failed to prepare indexing method %s\n", str);
+ free(str);
+ return NULL;
+ }
- }
+ if ( *m != INDEXING_NONE ) {
+ STATUS("Prepared indexing method %s\n", str);
+ }
+ free(str);
- if ( iprivs[n] == NULL ) return NULL;
+ if ( in != *m ) {
+ ERROR("Note: flags were altered to take into account "
+ "the information provided and/or the limitations "
+ "of the indexing method.\nPlease check the "
+ "methods listed above carefully.\n");
+ }
- str = indexer_str(indm[n]);
- STATUS("Prepared indexing method %i: %s\n", n, str);
- free(str);
+ return priv;
+}
+
+
+IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
+ struct detector *det, float *ltl,
+ int no_refine, const char *options)
+{
+ int i, n;
+ char **method_strings;
+ IndexingPrivate *ipriv;
+
+ /* Parse indexing methods */
+ n = assplode(method_list, ",", &method_strings, ASSPLODE_NONE);
+
+ IndexingMethod *methods = malloc(n * sizeof(IndexingMethod));
+ if ( methods == NULL ) {
+ ERROR("Failed to allocate indexing method list\n");
+ return NULL;
+ }
- if ( in != indm[n] ) {
- ERROR("Note: flags were altered to take into account "
- "the information provided and/or the limitations "
- "of the indexing method.\nPlease check the "
- "methods listed above carefully.\n");
+ for ( i=0; i<n; i++ ) {
+ methods[i] = get_indm_from_string(method_strings[i]);
+ if ( methods[i] == INDEXING_ERROR ) {
+ free(methods);
+ return NULL;
}
+ }
+
+ ipriv = malloc(sizeof(struct _indexingprivate));
+ if ( ipriv == NULL ) {
+ ERROR("Failed to allocate indexing data\n");
+ return NULL;
+ }
+
+ ipriv->engine_private = malloc((n+1) * sizeof(void *));
- for ( i=0; i<n; i++ ) {
- if ( indm[i] == indm[n] ) {
+ for ( i=0; i<n; i++ ) {
+
+ int j;
+
+ ipriv->engine_private[i] = prepare_method(&methods[i], cell,
+ det, ltl, options);
+ for ( j=0; j<i; j++ ) {
+ if ( methods[i] == methods[j] ) {
ERROR("Duplicate indexing method.\n");
ERROR("Have you specified some flags which "
"aren't accepted by one of your "
@@ -152,59 +204,66 @@ IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
}
}
+ }
+
+ ipriv->methods = methods;
+ ipriv->n_methods = n;
+ if ( cell != NULL ) {
+ ipriv->target_cell = cell_new_from_cell(cell);
+ } else {
+ ipriv->target_cell = NULL;
}
- iprivs[n] = NULL;
+ for ( i=0; i<4; i++ ) ipriv->tolerance[i] = ltl[i];
- return iprivs;
+ return ipriv;
}
-void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs)
+void cleanup_indexing(IndexingPrivate *ipriv)
{
- int n = 0;
+ int n;
- if ( indms == NULL ) return; /* Nothing to do */
- if ( privs == NULL ) return; /* Nothing to do */
+ if ( ipriv == NULL ) return; /* Nothing to do */
- while ( indms[n] != INDEXING_NONE ) {
+ for ( n=0; n<ipriv->n_methods; n++ ) {
- switch ( indms[n] & INDEXING_METHOD_MASK ) {
+ switch ( ipriv->methods[n] & INDEXING_METHOD_MASK ) {
case INDEXING_NONE :
break;
case INDEXING_DIRAX :
- dirax_cleanup(privs[n]);
+ dirax_cleanup(ipriv->engine_private[n]);
break;
case INDEXING_ASDF :
- asdf_cleanup(privs[n]);
+ asdf_cleanup(ipriv->engine_private[n]);
break;
case INDEXING_MOSFLM :
- mosflm_cleanup(privs[n]);
+ mosflm_cleanup(ipriv->engine_private[n]);
break;
case INDEXING_XDS :
- xds_cleanup(privs[n]);
+ xds_cleanup(ipriv->engine_private[n]);
break;
case INDEXING_FELIX :
- felix_cleanup(privs[n]);
+ felix_cleanup(ipriv->engine_private[n]);
break;
case INDEXING_DEBUG :
- free(privs[n]);
+ free(ipriv->engine_private[n]);
break;
case INDEXING_TAKETWO :
- taketwo_cleanup(privs[n]);
+ taketwo_cleanup(ipriv->engine_private[n]);
break;
default :
ERROR("Don't know how to clean up indexing method %i\n",
- indms[n]);
+ ipriv->methods[n]);
break;
}
@@ -213,8 +272,10 @@ void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs)
}
- free(indms);
- free(privs);
+ free(ipriv->methods);
+ free(ipriv->engine_private);
+ cell_free(ipriv->target_cell);
+ free(ipriv);
}
@@ -243,7 +304,7 @@ void map_all_peaks(struct image *image)
/* Return non-zero for "success" */
static int try_indexer(struct image *image, IndexingMethod indm,
- IndexingPrivate *ipriv)
+ IndexingPrivate *ipriv, void *mpriv)
{
int i, r;
int n_bad = 0;
@@ -254,19 +315,19 @@ static int try_indexer(struct image *image, IndexingMethod indm,
return 0;
case INDEXING_DIRAX :
- r = run_dirax(image, ipriv);
+ r = run_dirax(image, mpriv);
break;
case INDEXING_ASDF :
- r = run_asdf(image, ipriv);
+ r = run_asdf(image, mpriv);
break;
case INDEXING_MOSFLM :
- r = run_mosflm(image, ipriv);
+ r = run_mosflm(image, mpriv);
break;
case INDEXING_XDS :
- r = run_xds(image, ipriv);
+ r = run_xds(image, mpriv);
break;
case INDEXING_DEBUG :
@@ -274,11 +335,11 @@ static int try_indexer(struct image *image, IndexingMethod indm,
break;
case INDEXING_FELIX :
- r = felix_index(image, ipriv);
+ r = felix_index(image, mpriv);
break;
case INDEXING_TAKETWO :
- r = taketwo_index(image, ipriv);
+ r = taketwo_index(image, mpriv);
break;
default :
@@ -296,11 +357,35 @@ static int try_indexer(struct image *image, IndexingMethod indm,
crystal_set_mosaicity(cr, 0.0);
/* Prediction refinement if requested */
- if ( (indm & INDEXING_REFINE)
- && (refine_prediction(image, cr) != 0) )
- {
- crystal_set_user_flag(cr, 1);
- n_bad++;
+ if ( indm & INDEXING_REFINE ) {
+
+ UnitCell *out;
+
+ if ( refine_prediction(image, cr) ) {
+ crystal_set_user_flag(cr, 1);
+ n_bad++;
+ continue;
+ }
+
+ if ( (indm & INDEXING_CHECK_CELL_COMBINATIONS)
+ || (indm & INDEXING_CHECK_CELL_AXES) )
+ {
+
+ /* Check that the cell parameters are still
+ * within the tolerance */
+ out = match_cell(crystal_get_cell(cr),
+ ipriv->target_cell, 0,
+ ipriv->tolerance, 0);
+
+ if ( out == NULL ) {
+ crystal_set_user_flag(cr, 1);
+ n_bad++;
+ }
+
+ cell_free(out);
+
+ }
+
}
}
@@ -434,20 +519,18 @@ static int finished_retry(IndexingMethod indm, int r, struct image *image)
}
}
-void index_pattern(struct image *image,
- IndexingMethod *indms, IndexingPrivate **iprivs)
+void index_pattern(struct image *image, IndexingPrivate *ipriv)
{
- index_pattern_2(image, indms, iprivs, NULL);
+ index_pattern_2(image, ipriv, NULL);
}
-void index_pattern_2(struct image *image, IndexingMethod *indms,
- IndexingPrivate **iprivs, int *ping)
+void index_pattern_2(struct image *image, IndexingPrivate *ipriv, int *ping)
{
int n = 0;
ImageFeatureList *orig;
- if ( indms == NULL ) return;
+ if ( ipriv == NULL ) return;
map_all_peaks(image);
image->crystals = NULL;
@@ -455,7 +538,7 @@ void index_pattern_2(struct image *image, IndexingMethod *indms,
orig = image->features;
- while ( indms[n] != INDEXING_NONE ) {
+ for ( n=0; n<ipriv->n_methods; n++ ) {
int done = 0;
int r;
@@ -466,10 +549,11 @@ void index_pattern_2(struct image *image, IndexingMethod *indms,
do {
- r = try_indexer(image, indms[n], iprivs[n]);
+ r = try_indexer(image, ipriv->methods[n],
+ ipriv, ipriv->engine_private[n]);
success += r;
ntry++;
- done = finished_retry(indms[n], r, image);
+ done = finished_retry(ipriv->methods[n], r, image);
if ( ntry > 5 ) done = 1;
if ( ping != NULL ) (*ping)++;
@@ -481,11 +565,14 @@ void index_pattern_2(struct image *image, IndexingMethod *indms,
* crystals with a different indexing method) */
if ( success ) break;
- n++;
+ }
+ if ( n < ipriv->n_methods ) {
+ image->indexed_by = ipriv->methods[n];
+ } else {
+ image->indexed_by = INDEXING_NONE;
}
- image->indexed_by = indms[n];
image->features = orig;
}
@@ -656,100 +743,123 @@ char *indexer_str(IndexingMethod indm)
}
-IndexingMethod *build_indexer_list(const char *str)
+static IndexingMethod warn_method(const char *str)
+{
+ ERROR("Indexing method must contain exactly one engine name: '%s'\n",
+ str);
+ return INDEXING_ERROR;
+}
+
+
+IndexingMethod get_indm_from_string(const char *str)
{
int n, i;
- char **methods;
- IndexingMethod *list;
- int nmeth = 0;
+ char **bits;
+ IndexingMethod method = INDEXING_NONE;
+ int have_method = 0;
- n = assplode(str, ",-", &methods, ASSPLODE_NONE);
- list = malloc((n+1)*sizeof(IndexingMethod));
+ n = assplode(str, "-", &bits, ASSPLODE_NONE);
- nmeth = -1; /* So that the first method is #0 */
for ( i=0; i<n; i++ ) {
- if ( strcmp(methods[i], "dirax") == 0) {
- list[++nmeth] = INDEXING_DEFAULTS_DIRAX;
+ if ( strcmp(bits[i], "dirax") == 0) {
+ if ( have_method ) return warn_method(str);
+ method = INDEXING_DEFAULTS_DIRAX;
+ have_method = 1;
- } else if ( strcmp(methods[i], "asdf") == 0) {
- list[++nmeth] = INDEXING_DEFAULTS_ASDF;
+ } else if ( strcmp(bits[i], "asdf") == 0) {
+ if ( have_method ) return warn_method(str);
+ method = INDEXING_DEFAULTS_ASDF;
+ have_method = 1;
- } else if ( strcmp(methods[i], "mosflm") == 0) {
- list[++nmeth] = INDEXING_DEFAULTS_MOSFLM;
+ } else if ( strcmp(bits[i], "mosflm") == 0) {
+ if ( have_method ) return warn_method(str);
+ method = INDEXING_DEFAULTS_MOSFLM;
+ have_method = 1;
- } else if ( strcmp(methods[i], "xds") == 0) {
- list[++nmeth] = INDEXING_DEFAULTS_XDS;
+ } else if ( strcmp(bits[i], "xds") == 0) {
+ if ( have_method ) return warn_method(str);
+ method = INDEXING_DEFAULTS_XDS;
+ have_method = 1;
- } else if ( strcmp(methods[i], "felix") == 0) {
- list[++nmeth] = INDEXING_DEFAULTS_FELIX;
+ } else if ( strcmp(bits[i], "felix") == 0) {
+ if ( have_method ) return warn_method(str);
+ method = INDEXING_DEFAULTS_FELIX;
+ have_method = 1;
- } else if ( strcmp(methods[i], "taketwo") == 0) {
- list[++nmeth] = INDEXING_DEFAULTS_TAKETWO;
+ } else if ( strcmp(bits[i], "taketwo") == 0) {
+ if ( have_method ) return warn_method(str);
+ method = INDEXING_DEFAULTS_TAKETWO;
+ have_method = 1;
- } else if ( strcmp(methods[i], "none") == 0) {
- list[++nmeth] = INDEXING_NONE;
+ } else if ( strcmp(bits[i], "none") == 0) {
+ if ( have_method ) return warn_method(str);
+ method = INDEXING_NONE;
+ have_method = 1;
- } else if ( strcmp(methods[i], "simulation") == 0) {
- list[++nmeth] = INDEXING_SIMULATION;
- return list;
+ } else if ( strcmp(bits[i], "simulation") == 0) {
+ if ( have_method ) return warn_method(str);
+ method = INDEXING_SIMULATION;
+ return method;
- } else if ( strcmp(methods[i], "debug") == 0) {
- list[++nmeth] = INDEXING_DEBUG;
- return list;
+ } else if ( strcmp(bits[i], "debug") == 0) {
+ if ( have_method ) return warn_method(str);
+ method = INDEXING_DEBUG;
+ return method;
- } else if ( strcmp(methods[i], "raw") == 0) {
- list[nmeth] = set_raw(list[nmeth]);
+ } else if ( strcmp(bits[i], "raw") == 0) {
+ method = set_raw(method);
- } else if ( strcmp(methods[i], "bad") == 0) {
- list[nmeth] = set_bad(list[nmeth]);
+ } else if ( strcmp(bits[i], "bad") == 0) {
+ method = set_bad(method);
- } else if ( strcmp(methods[i], "comb") == 0) {
- list[nmeth] = set_comb(list[nmeth]); /* Default */
+ } else if ( strcmp(bits[i], "comb") == 0) {
+ method = set_comb(method); /* Default */
- } else if ( strcmp(methods[i], "axes") == 0) {
- list[nmeth] = set_axes(list[nmeth]);
+ } else if ( strcmp(bits[i], "axes") == 0) {
+ method = set_axes(method);
- } else if ( strcmp(methods[i], "latt") == 0) {
- list[nmeth] = set_lattice(list[nmeth]);
+ } else if ( strcmp(bits[i], "latt") == 0) {
+ method = set_lattice(method);
- } else if ( strcmp(methods[i], "nolatt") == 0) {
- list[nmeth] = set_nolattice(list[nmeth]);
+ } else if ( strcmp(bits[i], "nolatt") == 0) {
+ method = set_nolattice(method);
- } else if ( strcmp(methods[i], "cell") == 0) {
- list[nmeth] = set_cellparams(list[nmeth]);
+ } else if ( strcmp(bits[i], "cell") == 0) {
+ method = set_cellparams(method);
- } else if ( strcmp(methods[i], "nocell") == 0) {
- list[nmeth] = set_nocellparams(list[nmeth]);
+ } else if ( strcmp(bits[i], "nocell") == 0) {
+ method = set_nocellparams(method);
- } else if ( strcmp(methods[i], "retry") == 0) {
- list[nmeth] |= INDEXING_RETRY;
+ } else if ( strcmp(bits[i], "retry") == 0) {
+ method |= INDEXING_RETRY;
- } else if ( strcmp(methods[i], "noretry") == 0) {
- list[nmeth] &= ~INDEXING_RETRY;
+ } else if ( strcmp(bits[i], "noretry") == 0) {
+ method &= ~INDEXING_RETRY;
- } else if ( strcmp(methods[i], "multi") == 0) {
- list[nmeth] |= INDEXING_MULTI;
+ } else if ( strcmp(bits[i], "multi") == 0) {
+ method |= INDEXING_MULTI;
- } else if ( strcmp(methods[i], "nomulti") == 0) {
- list[nmeth] &= ~INDEXING_MULTI;
+ } else if ( strcmp(bits[i], "nomulti") == 0) {
+ method &= ~INDEXING_MULTI;
- } else if ( strcmp(methods[i], "refine") == 0) {
- list[nmeth] |= INDEXING_REFINE;
+ } else if ( strcmp(bits[i], "refine") == 0) {
+ method |= INDEXING_REFINE;
- } else if ( strcmp(methods[i], "norefine") == 0) {
- list[nmeth] &= ~INDEXING_REFINE;
+ } else if ( strcmp(bits[i], "norefine") == 0) {
+ method &= ~INDEXING_REFINE;
} else {
ERROR("Bad list of indexing methods: '%s'\n", str);
- return NULL;
+ return INDEXING_ERROR;
}
- free(methods[i]);
+ free(bits[i]);
}
- free(methods);
- list[++nmeth] = INDEXING_NONE;
+ free(bits);
+
+ if ( !have_method ) return warn_method(str);
- return list;
+ return method;
}
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index 33d568f6..2107ed80 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -81,6 +81,7 @@
* @INDEXING_DEBUG: Results injector for debugging
* @INDEXING_ASDF: Use in-built "asdf" indexer
* @INDEXING_TAKETWO: Use in-built "taketwo" indexer
+ * @INDEXING_ERROR: Special value for unrecognised indexing engine name
* @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
@@ -90,7 +91,7 @@
* @INDEXING_USE_LATTICE_TYPE: Use lattice type and centering information to
* guide the indexing process.
* @INDEXING_USE_CELL_PARAMETERS: Use the unit cell parameters to guide the
- * indexingprocess.
+ * indexing process.
* @INDEXING_RETRY: If the indexer doesn't succeed, delete the weakest peaks
* and try again.
* @INDEXING_MULTI: If the indexer succeeds, delete the peaks explained by the
@@ -115,6 +116,8 @@ typedef enum {
INDEXING_ASDF = 8,
INDEXING_TAKETWO = 9,
+ INDEXING_ERROR = 255, /* Unrecognised indexing engine */
+
/* Bits at the top of the IndexingMethod are flags which modify the
* behaviour of the indexer. */
INDEXING_CHECK_CELL_COMBINATIONS = 256,
@@ -143,28 +146,28 @@ extern "C" {
* IndexingPrivate:
*
* This is an opaque data structure containing information needed by the
- * indexing method.
+ * indexing system.
**/
-typedef void *IndexingPrivate;
+typedef struct _indexingprivate IndexingPrivate;
-extern IndexingMethod *build_indexer_list(const char *str);
+/* Convert indexing methods to/from text */
extern char *indexer_str(IndexingMethod indm);
+extern IndexingMethod get_indm_from_string(const char *method);
#include "detector.h"
#include "cell.h"
#include "image.h"
-extern IndexingPrivate **prepare_indexing(IndexingMethod *indm, UnitCell *cell,
- struct detector *det, float *ltl,
- const char *options);
+extern IndexingPrivate *setup_indexing(const char *methods, UnitCell *cell,
+ struct detector *det, float *ltl,
+ int no_refine, const char *options);
-extern void index_pattern(struct image *image,
- IndexingMethod *indms, IndexingPrivate **iprivs);
+extern void index_pattern(struct image *image, IndexingPrivate *ipriv);
-extern void index_pattern_2(struct image *image, IndexingMethod *indms,
- IndexingPrivate **iprivs, int *ping);
+extern void index_pattern_2(struct image *image, IndexingPrivate *ipriv,
+ int *ping);
-extern void cleanup_indexing(IndexingMethod *indms, IndexingPrivate **privs);
+extern void cleanup_indexing(IndexingPrivate *ipriv);
#ifdef __cplusplus
}
diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c
index 05f853eb..a8e6e119 100644
--- a/libcrystfel/src/mosflm.c
+++ b/libcrystfel/src/mosflm.c
@@ -725,7 +725,7 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm)
}
-int run_mosflm(struct image *image, IndexingPrivate *ipriv)
+int run_mosflm(struct image *image, void *ipriv)
{
struct mosflm_data *mosflm;
unsigned int opts;
@@ -843,8 +843,8 @@ int run_mosflm(struct image *image, IndexingPrivate *ipriv)
}
-IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell,
- struct detector *det, float *ltl)
+void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct detector *det, float *ltl)
{
struct mosflm_private *mp;
int need_cell = 0;
@@ -911,7 +911,7 @@ IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell,
}
-void mosflm_cleanup(IndexingPrivate *pp)
+void mosflm_cleanup(void *pp)
{
struct mosflm_private *p;
p = (struct mosflm_private *)pp;
diff --git a/libcrystfel/src/mosflm.h b/libcrystfel/src/mosflm.h
index 572741dd..39ec5390 100644
--- a/libcrystfel/src/mosflm.h
+++ b/libcrystfel/src/mosflm.h
@@ -3,13 +3,13 @@
*
* Invoke the DPS auto-indexing algorithm through MOSFLM
*
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
* Authors:
* 2010 Richard Kirian <rkirian@asu.edu>
- * 2012-2014 Thomas White <taw@physics.org>
+ * 2012-2014,2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -41,12 +41,12 @@
extern "C" {
#endif
-extern int run_mosflm(struct image *image, IndexingPrivate *ipriv);
+extern int run_mosflm(struct image *image, void *ipriv);
-extern IndexingPrivate *mosflm_prepare(IndexingMethod *indm, UnitCell *cell,
- struct detector *det, float *ltl);
+extern void *mosflm_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct detector *det, float *ltl);
-extern void mosflm_cleanup(IndexingPrivate *pp);
+extern void mosflm_cleanup(void *pp);
#ifdef __cplusplus
}
diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c
new file mode 100644
index 00000000..417465df
--- /dev/null
+++ b/libcrystfel/src/peakfinder8.c
@@ -0,0 +1,1195 @@
+/*
+ * peakfinder8.c
+ *
+ * The peakfinder8 algorithm
+ *
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2017 Valerio Mariani <valerio.mariani@desy.de>
+ * 2017 Anton Barty <anton.barty@desy.de>
+ * 2017 Oleksandr Yefanov <oleksandr.yefanov@desy.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <stdlib.h>
+
+#include "peakfinder8.h"
+
+
+// CrystFEL-only block 1
+struct radius_maps
+{
+ float **r_maps;
+ int n_rmaps;
+};
+
+
+struct peakfinder_mask
+{
+ char **masks;
+ int n_masks;
+};
+
+
+struct peakfinder_panel_data
+{
+ float **panel_data;
+ int *panel_h;
+ int *panel_w;
+ int num_panels;
+};
+// End of CrystFEL-only block 1
+
+
+struct radial_stats
+{
+ float *roffset;
+ float *rthreshold;
+ float *lthreshold;
+ float *rsigma;
+ int *rcount;
+ int n_rad_bins;
+};
+
+
+struct peakfinder_intern_data
+{
+ char *pix_in_peak_map;
+ int *infs;
+ int *inss;
+ int *peak_pixels;
+};
+
+
+struct peakfinder_peak_data
+{
+ int num_found_peaks;
+ int *npix;
+ float *com_fs;
+ float *com_ss;
+ int *com_index;
+ float *tot_i;
+ float *max_i;
+ float *sigma;
+ float *snr;
+};
+
+
+// CrystFEL-only block 2
+static struct radius_maps *compute_radius_maps(struct detector *det)
+{
+ int i, u, iss, ifs;
+ struct panel p;
+ struct radius_maps *rm = NULL;
+
+ rm = (struct radius_maps *)malloc(sizeof(struct radius_maps));
+ if ( rm == NULL ) {
+ return NULL;
+ }
+
+ rm->r_maps = (float **)malloc(det->n_panels*sizeof(float*));
+ if ( rm->r_maps == NULL ) {
+ free(rm);
+ return NULL;
+ }
+
+ rm->n_rmaps = det->n_panels;
+
+ for( i=0 ; i<det->n_panels ; i++ ) {
+
+ p = det->panels[i];
+ rm->r_maps[i] = (float *)malloc(p.h*p.w*sizeof(float));
+
+ if ( rm->r_maps[i] == NULL ) {
+ for ( u = 0; u<i; u++ ) {
+ free(rm->r_maps[u]);
+ }
+ free(rm);
+ return NULL;
+ }
+
+ for ( iss=0 ; iss<p.h ; iss++ ) {
+ for ( ifs=0; ifs<p.w; ifs++ ) {
+
+ int rmi;
+ int x,y;
+
+ rmi = ifs + p.w * iss;
+
+ x = (p.cnx + ifs * p.fsx + iss * p.ssx);
+ y = (p.cny + ifs * p.fsy + iss * p.ssy);
+
+ rm->r_maps[i][rmi] = sqrt(x * x + y * y);
+ }
+ }
+ }
+ return rm;
+}
+
+
+static void free_radius_maps(struct radius_maps *r_maps)
+{
+ int i;
+
+ for ( i=0 ; i<r_maps->n_rmaps ; i++ ) {
+ free(r_maps->r_maps[i]);
+ }
+ free(r_maps->r_maps);
+ free(r_maps);
+}
+
+
+static struct peakfinder_mask *create_peakfinder_mask(struct image *img,
+ struct radius_maps *rmps,
+ int min_res,
+ int max_res)
+{
+ int i;
+ struct peakfinder_mask *msk;
+
+ msk = (struct peakfinder_mask *)malloc(sizeof(struct peakfinder_mask));
+ msk->masks =(char **) malloc(img->det->n_panels*sizeof(char*));
+ msk->n_masks = img->det->n_panels;
+ for ( i=0; i<img->det->n_panels; i++) {
+
+ struct panel p;
+ int iss, ifs;
+
+ p = img->det->panels[i];
+
+ msk->masks[i] = (char *)calloc(p.w*p.h,sizeof(char));
+
+ for ( iss=0 ; iss<p.h ; iss++ ) {
+ for ( ifs=0 ; ifs<p.w ; ifs++ ) {
+
+ int idx;
+
+ idx = ifs + iss*p.w;
+
+ if ( rmps->r_maps[i][idx] < max_res
+ && rmps->r_maps[i][idx] > min_res ) {
+
+ if (! ( ( img->bad != NULL )
+ && ( img->bad[i] != NULL )
+ && ( img->bad[i][idx] != 0 ) ) ) {
+ msk->masks[i][idx] = 1;
+ }
+
+ }
+ }
+ }
+ }
+ return msk;
+}
+
+
+static void free_peakfinder_mask(struct peakfinder_mask * pfmask)
+{
+ int i;
+
+ for ( i=0 ; i<pfmask->n_masks ; i++ ) {
+ free(pfmask->masks[i]);
+ }
+ free(pfmask->masks);
+ free(pfmask);
+}
+
+
+static struct peakfinder_panel_data *allocate_panel_data(int num_panels)
+{
+
+ struct peakfinder_panel_data *pfdata;
+
+ pfdata = (struct peakfinder_panel_data *)malloc(sizeof(struct peakfinder_panel_data));
+ if ( pfdata == NULL ) {
+ return NULL;
+ }
+
+ pfdata->panel_h = (int *)malloc(num_panels*sizeof(int));
+ if ( pfdata->panel_h == NULL ) {
+ free(pfdata);
+ return NULL;
+ }
+
+ pfdata->panel_w = (int *)malloc(num_panels*sizeof(int));
+ if ( pfdata->panel_w == NULL ) {
+ free(pfdata->panel_h);
+ free(pfdata);
+ return NULL;
+ }
+
+ pfdata->panel_data = (float **)malloc(num_panels*sizeof(float*));
+ if ( pfdata->panel_data == NULL ) {
+ free(pfdata->panel_w);
+ free(pfdata->panel_h);
+ free(pfdata);
+ return NULL;
+ }
+
+ pfdata->num_panels = num_panels;
+
+ return pfdata;
+}
+
+
+static void free_panel_data(struct peakfinder_panel_data *pfdata)
+{
+ free(pfdata->panel_data);
+ free(pfdata->panel_w);
+ free(pfdata->panel_h);
+ free(pfdata);
+}
+
+
+static void compute_num_radial_bins(int w, int h, float *r_map, float *max_r)
+{
+ int ifs, iss;
+ int pidx;
+
+ for ( iss=0 ; iss<h ; iss++ ) {
+ for ( ifs=0 ; ifs<w ; ifs++ ) {
+ pidx = iss * w + ifs;
+ if ( r_map[pidx] > *max_r ) {
+ *max_r = r_map[pidx];
+ }
+ }
+ }
+}
+// End of CrystFEL-only block 2
+
+
+static struct radial_stats* allocate_radial_stats(int num_rad_bins)
+{
+ struct radial_stats* rstats;
+
+ rstats = (struct radial_stats *)malloc(sizeof(struct radial_stats));
+ if ( rstats == NULL ) {
+ return NULL;
+ }
+
+ rstats->roffset = (float *)malloc(num_rad_bins*sizeof(float));
+ if ( rstats->roffset == NULL ) {
+ free(rstats);
+ return NULL;
+ }
+
+ rstats->rthreshold = (float *)malloc(num_rad_bins*sizeof(float));
+ if ( rstats->rthreshold == NULL ) {
+ free(rstats);
+ free(rstats->roffset);
+ return NULL;
+ }
+
+ rstats->lthreshold = (float *)malloc(num_rad_bins*sizeof(float));
+ if ( rstats->lthreshold == NULL ) {
+ free(rstats);
+ free(rstats->rthreshold);
+ free(rstats->roffset);
+ return NULL;
+ }
+
+ rstats->rsigma = (float *)malloc(num_rad_bins*sizeof(float));
+ if ( rstats->rsigma == NULL ) {
+ free(rstats);
+ free(rstats->roffset);
+ free(rstats->rthreshold);
+ free(rstats->lthreshold);
+ return NULL;
+ }
+
+ rstats->rcount = (int *)malloc(num_rad_bins*sizeof(int));
+ if ( rstats->rcount == NULL ) {
+ free(rstats);
+ free(rstats->roffset);
+ free(rstats->rthreshold);
+ free(rstats->lthreshold);
+ free(rstats->rsigma);
+ return NULL;
+ }
+
+ rstats->n_rad_bins = num_rad_bins;
+
+ return rstats;
+}
+
+
+static void free_radial_stats(struct radial_stats *rstats)
+{
+ free(rstats->roffset);
+ free(rstats->rthreshold);
+ free(rstats->lthreshold);
+ free(rstats->rsigma);
+ free(rstats->rcount);
+ free(rstats);
+}
+
+
+static void fill_radial_bins(float *data,
+ int w,
+ int h,
+ float *r_map,
+ char *mask,
+ float *rthreshold,
+ float *lthreshold,
+ float *roffset,
+ float *rsigma,
+ int *rcount)
+{
+ int iss, ifs;
+ int pidx;
+
+ int curr_r;
+ float value;
+
+ for ( iss = 0; iss<h ; iss++ ) {
+ for ( ifs = 0; ifs<w ; ifs++ ) {
+ pidx = iss * w + ifs;
+ if ( mask[pidx] != 0 ) {
+ curr_r = (int)rint(r_map[pidx]);
+ value = data[pidx];
+ if ( value < rthreshold[curr_r ] && value>lthreshold[curr_r]) {
+ roffset[curr_r] += value;
+ rsigma[curr_r] += (value * value);
+ rcount[curr_r] += 1;
+ }
+ }
+ }
+ }
+}
+
+
+static void compute_radial_stats(float *rthreshold,
+ float *lthreshold,
+ float *roffset,
+ float *rsigma,
+ int *rcount,
+ int num_rad_bins,
+ float min_snr,
+ float acd_threshold)
+{
+ int ri;
+ float this_offset, this_sigma;
+
+ for ( ri=0 ; ri<num_rad_bins ; ri++ ) {
+
+ if ( rcount[ri] == 0 ) {
+ roffset[ri] = 0;
+ rsigma[ri] = 0;
+ rthreshold[ri] = 1e9;
+ lthreshold[ri] = -1e9;
+ } else {
+ this_offset = roffset[ri] / rcount[ri];
+ this_sigma = rsigma[ri] / rcount[ri] - (this_offset * this_offset);
+ if ( this_sigma >= 0 ) {
+ this_sigma = sqrt(this_sigma);
+ }
+
+ roffset[ri] = this_offset;
+ rsigma[ri] = this_sigma;
+ rthreshold[ri] = roffset[ri] + min_snr*rsigma[ri];
+ lthreshold[ri] = roffset[ri] - min_snr*rsigma[ri];
+
+ if ( rthreshold[ri] < acd_threshold ) {
+ rthreshold[ri] = acd_threshold;
+ }
+ }
+ }
+
+}
+
+
+struct peakfinder_peak_data *allocate_peak_data(int max_num_peaks)
+{
+ struct peakfinder_peak_data *pkdata;
+
+ pkdata = (struct peakfinder_peak_data*)malloc(sizeof(struct peakfinder_peak_data));
+ if ( pkdata == NULL ) {
+ return NULL;
+ }
+
+ pkdata->npix = (int *)malloc(max_num_peaks*sizeof(int));
+ if ( pkdata->npix == NULL ) {
+ free(pkdata);
+ free(pkdata->npix);
+ return NULL;
+ }
+
+ pkdata->com_fs = (float *)malloc(max_num_peaks*sizeof(float));
+ if ( pkdata->com_fs == NULL ) {
+ free(pkdata->npix);
+ free(pkdata);
+ return NULL;
+ }
+
+ pkdata->com_ss = (float *)malloc(max_num_peaks*sizeof(float));
+ if ( pkdata->com_ss == NULL ) {
+ free(pkdata->npix);
+ free(pkdata->com_fs);
+ free(pkdata);
+ return NULL;
+ }
+
+ pkdata->com_index = (int *)malloc(max_num_peaks*sizeof(int));
+ if ( pkdata->com_ss == NULL ) {
+ free(pkdata->npix);
+ free(pkdata->com_fs);
+ free(pkdata->com_ss);
+ free(pkdata);
+ return NULL;
+ }
+
+ pkdata->tot_i = (float *)malloc(max_num_peaks*sizeof(float));
+ if ( pkdata->tot_i == NULL ) {
+ free(pkdata->npix);
+ free(pkdata->com_fs);
+ free(pkdata->com_ss);
+ free(pkdata->com_index);
+ free(pkdata);
+ return NULL;
+ }
+
+ pkdata->max_i = (float *)malloc(max_num_peaks*sizeof(float));
+ if ( pkdata->max_i == NULL ) {
+ free(pkdata->npix);
+ free(pkdata->com_fs);
+ free(pkdata->com_ss);
+ free(pkdata->com_index);
+ free(pkdata->tot_i);
+ free(pkdata);
+ return NULL;
+ }
+
+ pkdata->sigma = (float *)malloc(max_num_peaks*sizeof(float));
+ if ( pkdata->sigma == NULL ) {
+ free(pkdata->npix);
+ free(pkdata->com_fs);
+ free(pkdata->com_ss);
+ free(pkdata->com_index);
+ free(pkdata->tot_i);
+ free(pkdata->max_i);
+ free(pkdata);
+ return NULL;
+ }
+
+ pkdata->snr = (float *)malloc(max_num_peaks*sizeof(float));
+ if ( pkdata->snr == NULL ) {
+ free(pkdata->npix);
+ free(pkdata->com_fs);
+ free(pkdata->com_ss);
+ free(pkdata->com_index);
+ free(pkdata->tot_i);
+ free(pkdata->max_i);
+ free(pkdata->sigma);
+ free(pkdata);
+ return NULL;
+ }
+
+ return pkdata;
+}
+
+
+static void free_peak_data(struct peakfinder_peak_data *pkdata) {
+ free(pkdata->npix);
+ free(pkdata->com_fs);
+ free(pkdata->com_ss);
+ free(pkdata->com_index);
+ free(pkdata->tot_i);
+ free(pkdata->max_i);
+ free(pkdata->sigma);
+ free(pkdata->snr);
+ free(pkdata);
+}
+
+
+static struct peakfinder_intern_data *allocate_peakfinder_intern_data(int data_size,
+ int max_pix_count)
+{
+
+ struct peakfinder_intern_data *intern_data;
+
+ intern_data = (struct peakfinder_intern_data *)malloc(sizeof(struct peakfinder_intern_data));
+ if ( intern_data == NULL ) {
+ return NULL;
+ }
+
+ intern_data->pix_in_peak_map =(char *)calloc(data_size, sizeof(char));
+ if ( intern_data->pix_in_peak_map == NULL ) {
+ free(intern_data);
+ return NULL;
+ }
+
+ intern_data->infs =(int *)calloc(data_size, sizeof(int));
+ if ( intern_data->infs == NULL ) {
+ free(intern_data->pix_in_peak_map);
+ free(intern_data);
+ return NULL;
+ }
+
+ intern_data->inss =(int *)calloc(data_size, sizeof(int));
+ if ( intern_data->inss == NULL ) {
+ free(intern_data->pix_in_peak_map);
+ free(intern_data->infs);
+ free(intern_data);
+ return NULL;
+ }
+
+ intern_data->peak_pixels =(int *)calloc(max_pix_count, sizeof(int));
+ if ( intern_data->peak_pixels == NULL ) {
+ free(intern_data->pix_in_peak_map);
+ free(intern_data->infs);
+ free(intern_data->inss);
+ free(intern_data);
+ return NULL;
+ }
+
+ return intern_data;
+}
+
+
+static void free_peakfinder_intern_data(struct peakfinder_intern_data *pfid)
+{
+ free(pfid->peak_pixels);
+ free(pfid->pix_in_peak_map);
+ free(pfid->infs);
+ free(pfid->inss);
+ free(pfid);
+}
+
+
+
+static void peak_search(int p,
+ struct peakfinder_intern_data *pfinter,
+ float *copy, char *mask, float *r_map,
+ float *rthreshold, float *roffset,
+ int *num_pix_in_peak, int asic_size_fs,
+ int asic_size_ss, int aifs, int aiss,
+ int num_pix_fs, float *sum_com_fs,
+ float *sum_com_ss, float *sum_i, int max_pix_count)
+{
+
+ int k, pi;
+ int curr_radius;
+ float curr_threshold;
+ int curr_fs;
+ int curr_ss;
+ float curr_i;
+
+ int search_fs[9] = { 0, -1, 0, 1, -1, 1, -1, 0, 1 };
+ int search_ss[9] = { 0, -1, -1, -1, 0, 0, 1, 1, 1 };
+ int search_n = 9;
+
+ // Loop through search pattern
+ for ( k=0; k<search_n; k++ ) {
+
+ if ( (pfinter->infs[p] + search_fs[k]) < 0 ) continue;
+ if ( (pfinter->infs[p] + search_fs[k]) >= asic_size_fs ) continue;
+ if ( (pfinter->inss[p] + search_ss[k]) < 0 ) continue;
+ if ( (pfinter->inss[p] + search_ss[k]) >= asic_size_ss ) continue;
+
+ // Neighbour point in big array
+ curr_fs = pfinter->infs[p] + search_fs[k] + aifs * asic_size_fs;
+ curr_ss = pfinter->inss[p] + search_ss[k] + aiss * asic_size_ss;
+ pi = curr_fs + curr_ss * num_pix_fs;
+
+ curr_radius = (int)rint(r_map[pi]);
+ curr_threshold = rthreshold[curr_radius];
+
+ // Above threshold?
+ if ( copy[pi] > curr_threshold
+ && pfinter->pix_in_peak_map[pi] == 0
+ && mask[pi] != 0 ) {
+
+ curr_i = copy[pi] - roffset[curr_radius];
+ *sum_i += curr_i;
+ *sum_com_fs += curr_i * ((float)curr_fs); // for center of mass x
+ *sum_com_ss += curr_i * ((float)curr_ss); // for center of mass y
+
+ pfinter->inss[*num_pix_in_peak] = pfinter->inss[p] + search_ss[k];
+ pfinter->infs[*num_pix_in_peak] = pfinter->infs[p] + search_fs[k];
+ pfinter->pix_in_peak_map[pi] = 1;
+ if ( *num_pix_in_peak < max_pix_count ) {
+ pfinter->peak_pixels[*num_pix_in_peak] = pi;
+ }
+ *num_pix_in_peak = *num_pix_in_peak + 1;
+ }
+ }
+}
+
+
+static void search_in_ring(int ring_width, int com_fs_int, int com_ss_int,
+ float *copy, float *r_map,
+ float *rthreshold, float *roffset,
+ char *pix_in_peak_map, char *mask, int asic_size_fs,
+ int asic_size_ss, int aifs, int aiss,
+ int num_pix_fs,float *local_sigma, float *local_offset,
+ float *background_max_i, int com_idx,
+ int local_bg_radius)
+{
+ int ssj, fsi;
+ float pix_radius;
+ int curr_fs, curr_ss;
+ int pi;
+ int curr_radius;
+ float curr_threshold;
+ float curr_i;
+
+ int np_sigma;
+ int np_counted;
+ int local_radius;
+
+ float sum_i;
+ float sum_i_squared;
+
+ ring_width = 2 * local_bg_radius;
+
+ sum_i = 0;
+ sum_i_squared = 0;
+ np_sigma = 0;
+ np_counted = 0;
+ local_radius = 0;
+
+ for ( ssj = -ring_width ; ssj<ring_width ; ssj++ ) {
+ for ( fsi = -ring_width ; fsi<ring_width ; fsi++ ) {
+
+ // Within-ASIC check
+ if ( (com_fs_int + fsi) < 0 ) continue;
+ if ( (com_fs_int + fsi) >= asic_size_fs ) continue;
+ if ( (com_ss_int + ssj) < 0 ) continue;
+ if ( (com_ss_int + ssj) >= asic_size_ss )
+ continue;
+
+ // Within outer ring check
+ pix_radius = sqrt(fsi * fsi + ssj * ssj);
+ if ( pix_radius>ring_width ) continue;
+
+ // Position of this point in data stream
+ curr_fs = com_fs_int + fsi + aifs * asic_size_fs;
+ curr_ss = com_ss_int + ssj + aiss * asic_size_ss;
+ pi = curr_fs + curr_ss * num_pix_fs;
+
+ curr_radius = (int)rint(r_map[pi]);
+ curr_threshold = rthreshold[curr_radius];
+
+ // Intensity above background ??? just intensity?
+ curr_i = copy[pi];
+
+ // Keep track of value and value-squared for offset and sigma calculation
+ if ( curr_i < curr_threshold && pix_in_peak_map[pi] == 0 && mask[pi] != 0 ) {
+
+ np_sigma++;
+ sum_i += curr_i;
+ sum_i_squared += (curr_i * curr_i);
+
+ if ( curr_i > *background_max_i ) {
+ *background_max_i = curr_i;
+ }
+ }
+ np_counted += 1;
+ }
+ }
+
+ // Calculate local background and standard deviation
+ if ( np_sigma != 0 ) {
+ *local_offset = sum_i / np_sigma;
+ *local_sigma = sum_i_squared / np_sigma - (*local_offset * *local_offset);
+ if (*local_sigma >= 0) {
+ *local_sigma = sqrt(*local_sigma);
+ } else {
+ *local_sigma = 0.01;
+ }
+ } else {
+ local_radius = (int)rint(r_map[(int)rint(com_idx)]);
+ *local_offset = roffset[local_radius];
+ *local_sigma = 0.01;
+ }
+}
+
+
+static void process_panel(int asic_size_fs, int asic_size_ss, int num_pix_fs,
+ int aiss, int aifs, float *rthreshold,
+ float *roffset, int *peak_count,
+ float *copy, struct peakfinder_intern_data *pfinter,
+ float *r_map, char *mask, int *npix, float *com_fs,
+ float *com_ss, int *com_index, float *tot_i,
+ float *max_i, float *sigma, float *snr,
+ int min_pix_count, int max_pix_count,
+ int local_bg_radius, float min_snr, int max_n_peaks)
+{
+ int pxss, pxfs;
+ int num_pix_in_peak;
+
+ // Loop over pixels within a module
+ for ( pxss=1 ; pxss<asic_size_ss-1 ; pxss++ ) {
+ for ( pxfs=1 ; pxfs<asic_size_fs-1 ; pxfs++ ) {
+
+ float curr_thresh;
+ int pxidx;
+ int curr_rad;
+
+ pxidx = (pxss + aiss * asic_size_ss) * num_pix_fs +
+ pxfs + aifs * asic_size_fs;
+
+ curr_rad = (int)rint(r_map[pxidx]);
+ curr_thresh = rthreshold[curr_rad];
+
+ if ( copy[pxidx] > curr_thresh
+ && pfinter->pix_in_peak_map[pxidx] == 0
+ && mask[pxidx] != 0 ) { //??? not sure if needed
+
+ // This might be the start of a new peak - start searching
+ float sum_com_fs, sum_com_ss;
+ float sum_i;
+ float peak_com_fs, peak_com_ss;
+ float peak_com_fs_int, peak_com_ss_int;
+ float peak_tot_i, pk_tot_i_raw;
+ float peak_max_i, pk_max_i_raw;
+ float peak_snr;
+ float local_sigma, local_offset;
+ float background_max_i;
+ int lt_num_pix_in_pk;
+ int ring_width;
+ int peak_idx;
+ int com_idx;
+ int p;
+
+ pfinter->infs[0] = pxfs;
+ pfinter->inss[0] = pxss;
+ pfinter->peak_pixels[0] = pxidx;
+ num_pix_in_peak = 0; //y 1;
+
+ sum_i = 0;
+ sum_com_fs = 0;
+ sum_com_ss = 0;
+
+ // Keep looping until the pixel count within this peak does not change
+ do {
+ lt_num_pix_in_pk = num_pix_in_peak;
+
+ // Loop through points known to be within this peak
+ for ( p=0; p<=num_pix_in_peak; p++ ) { //changed from 1 to 0 by O.Y.
+ peak_search(p,
+ pfinter, copy, mask,
+ r_map,
+ rthreshold,
+ roffset,
+ &num_pix_in_peak,
+ asic_size_fs,
+ asic_size_ss,
+ aifs, aiss,
+ num_pix_fs,
+ &sum_com_fs,
+ &sum_com_ss,
+ &sum_i,
+ max_pix_count);
+ }
+
+ } while ( lt_num_pix_in_pk != num_pix_in_peak );
+
+ // Too many or too few pixels means ignore this 'peak'; move on now
+ if ( num_pix_in_peak < min_pix_count || num_pix_in_peak > max_pix_count ) continue;
+
+ // If for some reason sum_i is 0 - it's better to skip
+ if ( fabs(sum_i) < 1e-10 ) continue;
+
+ // Calculate center of mass for this peak from initial peak search
+ peak_com_fs = sum_com_fs / fabs(sum_i);
+ peak_com_ss = sum_com_ss / fabs(sum_i);
+
+ com_idx = (int)rint(peak_com_fs) + (int)rint(peak_com_ss) * num_pix_fs;
+
+ peak_com_fs_int = (int)rint(peak_com_fs) - aifs * asic_size_fs;
+ peak_com_ss_int = (int)rint(peak_com_ss) - aiss * asic_size_ss;
+
+ // Calculate the local signal-to-noise ratio and local background in an annulus around
+ // this peak (excluding pixels which look like they might be part of another peak)
+ local_sigma = 0.0;
+ local_offset = 0.0;
+ background_max_i = 0.0;
+
+ ring_width = 2 * local_bg_radius;
+
+ search_in_ring(ring_width, peak_com_fs_int,
+ peak_com_ss_int,
+ copy, r_map, rthreshold,
+ roffset,
+ pfinter->pix_in_peak_map,
+ mask, asic_size_fs,
+ asic_size_ss,
+ aifs, aiss,
+ num_pix_fs,
+ &local_sigma,
+ &local_offset,
+ &background_max_i,
+ com_idx, local_bg_radius);
+
+ // Re-integrate (and re-centroid) peak using local background estimates
+ peak_tot_i = 0;
+ pk_tot_i_raw = 0;
+ peak_max_i = 0;
+ pk_max_i_raw = 0;
+ sum_com_fs = 0;
+ sum_com_ss = 0;
+
+ for ( peak_idx = 0 ;
+ peak_idx < num_pix_in_peak && peak_idx < max_pix_count ;
+ peak_idx++ ) {
+
+ int curr_idx;
+ float curr_i;
+ float curr_i_raw;
+ int curr_fs, curr_ss;
+
+ curr_idx = pfinter->peak_pixels[peak_idx];
+ curr_i_raw = copy[curr_idx];
+ curr_i = curr_i_raw - local_offset;
+ peak_tot_i += curr_i;
+ pk_tot_i_raw += curr_i_raw;
+
+ // Remember that curr_idx = curr_fs + curr_ss*num_pix_fs
+ curr_fs = curr_idx % num_pix_fs;
+ curr_ss = curr_idx / num_pix_fs;
+ sum_com_fs += curr_i * ((float)curr_fs);
+ sum_com_ss += curr_i * ((float)curr_ss);
+
+ if ( curr_i_raw > pk_max_i_raw ) pk_max_i_raw = curr_i_raw;
+ if ( curr_i > peak_max_i ) peak_max_i = curr_i;
+ }
+
+
+ // This CAN happen! Better to skip...
+ if ( fabs(peak_tot_i) < 1e-10 ) continue;
+
+ peak_com_fs = sum_com_fs / fabs(peak_tot_i);
+ peak_com_ss = sum_com_ss / fabs(peak_tot_i);
+
+ // Calculate signal-to-noise and apply SNR criteria
+ if ( fabs(local_sigma) > 1e-10 ) {
+ peak_snr = peak_tot_i / local_sigma;
+ } else {
+ peak_snr = 0;
+ }
+
+ if (peak_snr < min_snr) continue;
+
+ // Is the maximum intensity in the peak enough above intensity in background region to
+ // be a peak and not noise? The more pixels there are in the peak, the more relaxed we
+ // are about this criterion
+ //f_background_thresh = background_max_i - local_offset; //!!! Ofiget'! If I uncomment
+ // if (peak_max_i < f_background_thresh) { // these lines the result is
+ // different!
+ if (peak_max_i < background_max_i - local_offset) continue;
+
+ // This is a peak? If so, add info to peak list
+ if ( num_pix_in_peak >= min_pix_count
+ && num_pix_in_peak <= max_pix_count ) {
+
+ // Bragg peaks in the mask
+ for ( peak_idx = 0 ;
+ peak_idx < num_pix_in_peak &&
+ peak_idx < max_pix_count ;
+ peak_idx++ ) {
+ pfinter->pix_in_peak_map[pfinter->peak_pixels[peak_idx]] = 2;
+ }
+
+ int peak_com_idx;
+ peak_com_idx = (int)rint(peak_com_fs) + (int)rint(peak_com_ss) *
+ num_pix_fs;
+ // Remember peak information
+ if ( *peak_count < max_n_peaks ) {
+
+ int pidx;
+ pidx = *peak_count;
+
+ npix[pidx] = num_pix_in_peak;
+ com_fs[pidx] = peak_com_fs;
+ com_ss[pidx] = peak_com_ss;
+ com_index[pidx] = peak_com_idx;
+ tot_i[pidx] = peak_tot_i;
+ max_i[pidx] = peak_max_i;
+ sigma[pidx] = local_sigma;
+ snr[pidx] = peak_snr;
+ }
+ *peak_count += 1;
+ }
+ }
+ }
+ }
+}
+
+
+static int peakfinder8_base(float *roffset, float *rthreshold,
+ float *data, char *mask, float *r_map,
+ int asic_size_fs, int num_asics_fs,
+ int asic_size_ss, int num_asics_ss,
+ int max_n_peaks, int *num_found_peaks,
+ int *npix, float *com_fs,
+ float *com_ss, int *com_index, float *tot_i,
+ float *max_i, float *sigma, float *snr,
+ int min_pix_count, int max_pix_count,
+ int local_bg_radius, float min_snr,
+ char* outliersMask)
+{
+
+ int num_pix_fs, num_pix_ss, num_pix_tot;
+ int aifs, aiss;
+ int peak_count;
+ struct peakfinder_intern_data *pfinter;
+
+ num_pix_fs = asic_size_fs * num_asics_fs;
+ num_pix_ss = asic_size_ss * num_asics_ss;
+ num_pix_tot = num_pix_fs * num_pix_ss;
+
+ pfinter = allocate_peakfinder_intern_data(num_pix_tot, max_pix_count);
+ if ( pfinter == NULL ) {
+ return 1;
+ }
+
+ peak_count = 0;
+
+ // Loop over modules (nxn array)
+ for ( aiss=0 ; aiss<num_asics_ss ; aiss++ ) {
+ for ( aifs=0 ; aifs<num_asics_fs ; aifs++ ) { // ??? to change to proper panels need
+ process_panel(asic_size_fs, asic_size_ss, num_pix_fs, // change copy, mask, r_map
+ aiss, aifs, rthreshold, roffset,
+ &peak_count, data, pfinter, r_map, mask,
+ npix, com_fs, com_ss, com_index, tot_i,
+ max_i, sigma, snr, min_pix_count,
+ max_pix_count, local_bg_radius, min_snr,
+ max_n_peaks);
+ }
+ }
+ *num_found_peaks = peak_count;
+
+ if (outliersMask != NULL) {
+ memcpy(outliersMask, pfinter->pix_in_peak_map, num_pix_tot*sizeof(char));
+ }
+
+ free_peakfinder_intern_data(pfinter);
+
+ return 0;
+}
+
+
+int peakfinder8(struct image *img, int max_n_peaks,
+ float threshold, float min_snr,
+ int min_pix_count, int max_pix_count,
+ int local_bg_radius, int min_res,
+ int max_res, int use_saturated)
+{
+ struct radius_maps *rmaps;
+ struct peakfinder_mask *pfmask;
+ struct peakfinder_panel_data *pfdata;
+ struct radial_stats *rstats;
+ struct peakfinder_peak_data *pkdata;
+ int num_rad_bins;
+ int pi;
+ int i, it_counter;
+ int num_found_peaks;
+ int remaining_max_num_peaks;
+ int iterations;
+ float max_r;
+
+ iterations = 5;
+
+ if ( img-> det == NULL) {
+ return 1;
+ }
+
+ rmaps = compute_radius_maps(img->det);
+ if ( rmaps == NULL ) {
+ return 1;
+ }
+
+ pfmask = create_peakfinder_mask(img, rmaps, min_res, max_res);
+ if ( pfmask == NULL ) {
+ free_radius_maps(rmaps);
+ return 1;
+ }
+
+ pfdata = allocate_panel_data(img->det->n_panels);
+ if ( pfdata == NULL) {
+ free_radius_maps(rmaps);
+ free_peakfinder_mask(pfmask);
+ return 1;
+ }
+
+ for ( pi=0 ; pi<img->det->n_panels ; pi++ ) {
+ pfdata->panel_h[pi] = img->det->panels[pi].h;
+ pfdata->panel_w[pi] = img->det->panels[pi].w;
+ pfdata->panel_data[pi] = img->dp[pi];
+ pfdata->num_panels = img->det->n_panels;
+ }
+
+ max_r = -1e9;
+
+ for ( pi=0 ; pi<pfdata->num_panels ; pi++ ) {
+
+ compute_num_radial_bins(pfdata->panel_w[pi],
+ pfdata->panel_h[pi],
+ rmaps->r_maps[pi],
+ &max_r);
+ }
+
+ num_rad_bins = (int)ceil(max_r) + 1;
+
+ rstats = allocate_radial_stats(num_rad_bins);
+ if ( rstats == NULL ) {
+ free_radius_maps(rmaps);
+ free_peakfinder_mask(pfmask);
+ free_panel_data(pfdata);
+ return 1;
+ }
+
+ for ( i=0 ; i<rstats->n_rad_bins ; i++) {
+ rstats->rthreshold[i] = 1e9;
+ }
+
+ for ( it_counter=0 ; it_counter<iterations ; it_counter++ ) {
+
+ for ( i=0; i<num_rad_bins; i++ ) {
+ rstats->roffset[i] = 0;
+ rstats->rsigma[i] = 0;
+ rstats->rcount[i] = 0;
+ }
+
+ for ( pi=0 ; pi<pfdata->num_panels ; pi++ ) {
+
+ fill_radial_bins(pfdata->panel_data[pi],
+ pfdata->panel_w[pi],
+ pfdata->panel_h[pi],
+ rmaps->r_maps[pi],
+ pfmask->masks[pi],
+ rstats->rthreshold,
+ rstats->lthreshold,
+ rstats->roffset,
+ rstats->rsigma,
+ rstats->rcount);
+
+ }
+
+ compute_radial_stats(rstats->rthreshold,
+ rstats->lthreshold,
+ rstats->roffset,
+ rstats->rsigma,
+ rstats->rcount,
+ num_rad_bins,
+ min_snr,
+ threshold);
+
+ }
+
+ pkdata = allocate_peak_data(max_n_peaks);
+ if ( pkdata == NULL ) {
+ free_radius_maps(rmaps);
+ free_peakfinder_mask(pfmask);
+ free_panel_data(pfdata);
+ free_radial_stats(rstats);
+ return 1;
+ }
+
+ remaining_max_num_peaks = max_n_peaks;
+
+ for ( pi=0 ; pi<img->det->n_panels ; pi++) {
+
+ int peaks_to_add;
+ int pki;
+ int ret;
+
+ num_found_peaks = 0;
+
+ if ( img->det->panels[pi].no_index ) {
+ continue;
+ }
+
+ ret = peakfinder8_base(rstats->roffset,
+ rstats->rthreshold,
+ pfdata->panel_data[pi],
+ pfmask->masks[pi],
+ rmaps->r_maps[pi],
+ pfdata->panel_w[pi], 1,
+ pfdata->panel_h[pi], 1,
+ max_n_peaks,
+ &num_found_peaks,
+ pkdata->npix,
+ pkdata->com_fs,
+ pkdata->com_ss,
+ pkdata->com_index,
+ pkdata->tot_i,
+ pkdata->max_i,
+ pkdata->sigma,
+ pkdata->snr,
+ min_pix_count,
+ max_pix_count,
+ local_bg_radius,
+ min_snr,
+ NULL);
+
+ if ( ret != 0 ) {
+ free_radius_maps(rmaps);
+ free_peakfinder_mask(pfmask);
+ free_panel_data(pfdata);
+ free_radial_stats(rstats);
+ return 1;
+ }
+
+ peaks_to_add = num_found_peaks;
+
+ if ( num_found_peaks > remaining_max_num_peaks ) {
+ peaks_to_add = remaining_max_num_peaks;
+ }
+
+ remaining_max_num_peaks -= peaks_to_add;
+
+ for ( pki=0 ; pki<peaks_to_add ; pki++ ) {
+
+ struct panel *p;
+
+ p = &img->det->panels[pi];
+
+ img->num_peaks += 1;
+ if ( pkdata->max_i[pki] > p->max_adu ) {
+ img->num_saturated_peaks++;
+ if ( !use_saturated ) {
+ continue;
+ }
+ }
+
+ image_add_feature(img->features,
+ pkdata->com_fs[pki]+0.5,
+ pkdata->com_ss[pki]+0.5,
+ p,
+ img,
+ pkdata->tot_i[pki],
+ NULL);
+ }
+ }
+
+ free_radius_maps(rmaps);
+ free_peakfinder_mask(pfmask);
+ free_panel_data(pfdata);
+ free_radial_stats(rstats);
+ free_peak_data(pkdata);
+ return 0;
+}
diff --git a/libcrystfel/src/peakfinder8.h b/libcrystfel/src/peakfinder8.h
new file mode 100644
index 00000000..483eebdf
--- /dev/null
+++ b/libcrystfel/src/peakfinder8.h
@@ -0,0 +1,50 @@
+/*
+ * peakfinder8.h
+ *
+ * The peakfinder8 algorithm
+ *
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2017 Valerio Mariani <valerio.mariani@desy.de>
+ * 2017 Anton Barty <anton.barty@desy.de>
+ * 2017 Oleksandr Yefanov <oleksandr.yefanov@desy.de>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef PEAKFINDER8_H
+#define PEAKFINDER8_H
+
+#include "image.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern int peakfinder8(struct image *img, int max_n_peaks,
+ float threshold, float min_snr,
+ int mix_pix_count, int max_pix_count,
+ int local_bg_radius, int min_res,
+ int max_res, int use_saturated);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // PEAKFINDER8_H
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index cb2a9f5a..28c79538 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -3,7 +3,7 @@
*
* Peak search and other image analysis
*
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
@@ -12,6 +12,7 @@
* 2012 Kenneth Beyerlein <kenneth.beyerlein@desy.de>
* 2011 Andrew Martin <andrew.martin@desy.de>
* 2011 Richard Kirian
+ * 2017 Valerio Mariani <valerio.mariani@desy.de>
*
* This file is part of CrystFEL.
*
@@ -52,6 +53,7 @@
#include "reflist-utils.h"
#include "cell-utils.h"
#include "geometry.h"
+#include "peakfinder8.h"
static int cull_peaks_in_panel(struct image *image, struct panel *p)
@@ -517,8 +519,8 @@ static void search_peaks_in_panel(struct image *image, float threshold,
void search_peaks(struct image *image, float threshold, float min_gradient,
- float min_snr, double ir_inn, double ir_mid, double ir_out,
- int use_saturated)
+ float min_snr, double ir_inn, double ir_mid,
+ double ir_out, int use_saturated)
{
int i;
@@ -541,6 +543,26 @@ void search_peaks(struct image *image, float threshold, float min_gradient,
}
+int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
+ float threshold, float min_snr,
+ int min_pix_count, int max_pix_count,
+ int local_bg_radius, int min_res,
+ int max_res, int use_saturated)
+{
+ if ( image->features != NULL ) {
+ image_feature_list_free(image->features);
+ }
+ image->features = image_feature_list_new();
+ image->num_peaks = 0;
+ image->num_saturated_peaks = 0;
+
+ return peakfinder8(image, max_n_peaks, threshold, min_snr,
+ min_pix_count, max_pix_count,
+ local_bg_radius, min_res,
+ max_res, use_saturated);
+}
+
+
int peak_sanity_check(struct image *image, Crystal **crystals, int n_cryst)
{
int n_feat = 0;
diff --git a/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
index ba724419..a5095127 100644
--- a/libcrystfel/src/peaks.h
+++ b/libcrystfel/src/peaks.h
@@ -3,11 +3,12 @@
*
* Peak search and other image analysis
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
* 2010-2015 Thomas White <taw@physics.org>
+ * 2017 Valerio Mariani <valerio.mariani@desy.de>
*
* This file is part of CrystFEL.
*
@@ -45,9 +46,14 @@ extern "C" {
extern int *make_BgMask(struct image *image, struct panel *p, double ir_inn);
extern void search_peaks(struct image *image, float threshold,
- float min_gradient, float min_snr,
- double ir_inn, double ir_mid, double ir_out,
- int use_saturated);
+ float min_gradient, float min_snr, double ir_inn,
+ double ir_mid, double ir_out, int use_saturated);
+
+extern int search_peaks_peakfinder8(struct image *image, int max_n_peaks,
+ float threshold, float min_snr,
+ int mix_pix_count, int max_pix_count,
+ int local_bg_radius, int min_res,
+ int max_res, int use_saturated);
extern int peak_sanity_check(struct image *image, Crystal **crystals,
int n_cryst);
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
index 4621f4f4..380a1b8a 100644
--- a/libcrystfel/src/reflist-utils.c
+++ b/libcrystfel/src/reflist-utils.c
@@ -27,8 +27,6 @@
*
*/
-#define _ISOC99_SOURCE
-#define _GNU_SOURCE
#include <math.h>
#include <stdio.h>
#include <assert.h>
diff --git a/libcrystfel/src/statistics.c b/libcrystfel/src/statistics.c
index 56273fdb..ccf35194 100644
--- a/libcrystfel/src/statistics.c
+++ b/libcrystfel/src/statistics.c
@@ -30,7 +30,6 @@
#include <config.h>
#endif
-#define _ISOC99_SOURCE
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index e858172e..fb4b0c70 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -806,8 +806,8 @@ static int write_crystal(Stream *st, Crystal *cr, int include_reflections)
}
-int write_chunk(Stream *st, struct image *i, struct hdfile *hdfile,
- int include_peaks, int include_reflections, struct event* ev)
+int write_chunk(Stream *st, struct image *i, struct imagefile *imfile,
+ int include_peaks, int include_reflections, struct event *ev)
{
int j;
char *indexer;
@@ -832,7 +832,7 @@ int write_chunk(Stream *st, struct image *i, struct hdfile *hdfile,
fprintf(st->fh, "beam_divergence = %.2e rad\n", i->div);
fprintf(st->fh, "beam_bandwidth = %.2e (fraction)\n", i->bw);
- copy_hdf5_fields(hdfile, i->copyme, st->fh, ev);
+ imagefile_copy_fields(imfile, i->copyme, st->fh, ev);
if ( i->det != NULL ) {
@@ -1196,13 +1196,9 @@ int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf)
}
if ( strncmp(line, "indexed_by = ", 13) == 0 ) {
- IndexingMethod *list;
- list = build_indexer_list(line+13);
- if ( list == NULL ) {
+ image->indexed_by = get_indm_from_string(line+13);
+ if ( image->indexed_by == INDEXING_ERROR ) {
ERROR("Failed to read indexer list\n");
- } else {
- image->indexed_by = list[0];
- free(list);
}
}
diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h
index ff8628c0..764e3e36 100644
--- a/libcrystfel/src/stream.h
+++ b/libcrystfel/src/stream.h
@@ -3,11 +3,11 @@
*
* Stream tools
*
- * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2014 Thomas White <taw@physics.org>
+ * 2010-2017 Thomas White <taw@physics.org>
* 2014 Valerio Mariani
* 2011 Andrew Aquila
*
@@ -39,6 +39,7 @@
struct image;
struct hdfile;
struct event;
+struct imagefile;
#include "cell.h"
#define GEOM_START_MARKER "----- Begin geometry file -----"
@@ -106,10 +107,15 @@ extern int read_chunk(Stream *st, struct image *image);
extern int read_chunk_2(Stream *st, struct image *image,
StreamReadFlags srf);
-extern int write_chunk(Stream *st, struct image *image, struct hdfile *hdfile,
+extern int write_chunk(Stream *st, struct image *image, struct imagefile *imfile,
int include_peaks, int include_reflections,
struct event *ev);
+extern int write_chunk_2(Stream *st, struct image *image,
+ struct imagefile *imfile,
+ int include_peaks, int include_reflections,
+ struct event *ev);
+
extern void write_command(Stream *st, int argc, char *argv[]);
extern void write_geometry_file(Stream *st, const char *geom_filename);
extern int rewind_stream(Stream *st);
diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c
index 461da0cd..90d81a71 100644
--- a/libcrystfel/src/utils.c
+++ b/libcrystfel/src/utils.c
@@ -387,144 +387,6 @@ int poisson_noise(gsl_rng *rng, double expected)
}
-/**
- * SECTION:quaternion
- * @short_description: Simple quaternion handling
- * @title: Quaternion
- * @section_id:
- * @see_also:
- * @include: "utils.h"
- * @Image:
- *
- * There is a simple quaternion structure in CrystFEL. At the moment, it is
- * only used when simulating patterns, as an argument to cell_rotate() to
- * orient the unit cell.
- */
-
-/**
- * quaternion_modulus:
- * @q: A %quaternion
- *
- * If a quaternion represents a pure rotation, its modulus should be unity.
- *
- * Returns: the modulus of the given quaternion.
- **/
-double quaternion_modulus(struct quaternion q)
-{
- return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z);
-}
-
-
-/**
- * normalise_quaternion:
- * @q: A %quaternion
- *
- * Rescales the quaternion such that its modulus is unity.
- *
- * Returns: the normalised version of @q
- **/
-struct quaternion normalise_quaternion(struct quaternion q)
-{
- double mod;
- struct quaternion r;
-
- mod = quaternion_modulus(q);
-
- r.w = q.w / mod;
- r.x = q.x / mod;
- r.y = q.y / mod;
- r.z = q.z / mod;
-
- return r;
-}
-
-
-/**
- * random_quaternion:
- * @rng: A GSL random number generator to use
- *
- * Returns: a randomly generated, normalised, quaternion.
- **/
-struct quaternion random_quaternion(gsl_rng *rng)
-{
- struct quaternion q;
-
- q.w = 2.0*gsl_rng_uniform(rng) - 1.0;
- q.x = 2.0*gsl_rng_uniform(rng) - 1.0;
- q.y = 2.0*gsl_rng_uniform(rng) - 1.0;
- q.z = 2.0*gsl_rng_uniform(rng) - 1.0;
- q = normalise_quaternion(q);
-
- return q;
-}
-
-
-/**
- * quaternion_valid:
- * @q: A %quaternion
- *
- * Checks if the given quaternion is normalised.
- *
- * This function performs a nasty floating point comparison of the form
- * <code>(modulus > 0.999) && (modulus < 1.001)</code>, and so should not be
- * relied upon to spot anything other than the most obvious input error.
- *
- * Returns: 1 if the quaternion is normalised, 0 if not.
- **/
-int quaternion_valid(struct quaternion q)
-{
- double qmod;
-
- qmod = quaternion_modulus(q);
-
- /* Modulus = 1 to within some tolerance?
- * Nasty allowance for floating-point accuracy follows... */
- if ( (qmod > 0.999) && (qmod < 1.001) ) return 1;
-
- return 0;
-}
-
-
-/**
- * quat_rot
- * @q: A vector (in the form of a "struct rvec")
- * @z: A %quaternion
- *
- * Rotates a vector according to a quaternion.
- *
- * Returns: A rotated version of @p.
- **/
-struct rvec quat_rot(struct rvec q, struct quaternion z)
-{
- struct rvec res;
- double t01, t02, t03, t11, t12, t13, t22, t23, t33;
-
- t01 = z.w*z.x;
- t02 = z.w*z.y;
- t03 = z.w*z.z;
- t11 = z.x*z.x;
- t12 = z.x*z.y;
- t13 = z.x*z.z;
- t22 = z.y*z.y;
- t23 = z.y*z.z;
- t33 = z.z*z.z;
-
- res.u = (1.0 - 2.0 * (t22 + t33)) * q.u
- + (2.0 * (t12 + t03)) * q.v
- + (2.0 * (t13 - t02)) * q.w;
-
- res.v = (2.0 * (t12 - t03)) * q.u
- + (1.0 - 2.0 * (t11 + t33)) * q.v
- + (2.0 * (t01 + t23)) * q.w;
-
- res.w = (2.0 * (t02 + t13)) * q.u
- + (2.0 * (t23 - t01)) * q.v
- + (1.0 - 2.0 * (t11 + t22)) * q.w;
-
- return res;
-}
-
-
/* Return non-zero if c is in delims */
static int assplode_isdelim(const char c, const char *delims)
{
@@ -691,3 +553,141 @@ void utils_fudge_gslcblas()
{
STATUS("%p\n", cblas_sgemm);
}
+
+
+/**
+ * SECTION:quaternion
+ * @short_description: Simple quaternion handling
+ * @title: Quaternion
+ * @section_id:
+ * @see_also:
+ * @include: "utils.h"
+ * @Image:
+ *
+ * There is a simple quaternion structure in CrystFEL. At the moment, it is
+ * only used when simulating patterns, as an argument to cell_rotate() to
+ * orient the unit cell.
+ */
+
+/**
+ * quaternion_modulus:
+ * @q: A %quaternion
+ *
+ * If a quaternion represents a pure rotation, its modulus should be unity.
+ *
+ * Returns: the modulus of the given quaternion.
+ **/
+double quaternion_modulus(struct quaternion q)
+{
+ return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z);
+}
+
+
+/**
+ * normalise_quaternion:
+ * @q: A %quaternion
+ *
+ * Rescales the quaternion such that its modulus is unity.
+ *
+ * Returns: the normalised version of @q
+ **/
+struct quaternion normalise_quaternion(struct quaternion q)
+{
+ double mod;
+ struct quaternion r;
+
+ mod = quaternion_modulus(q);
+
+ r.w = q.w / mod;
+ r.x = q.x / mod;
+ r.y = q.y / mod;
+ r.z = q.z / mod;
+
+ return r;
+}
+
+
+/**
+ * random_quaternion:
+ * @rng: A GSL random number generator to use
+ *
+ * Returns: a randomly generated, normalised, quaternion.
+ **/
+struct quaternion random_quaternion(gsl_rng *rng)
+{
+ struct quaternion q;
+
+ q.w = 2.0*gsl_rng_uniform(rng) - 1.0;
+ q.x = 2.0*gsl_rng_uniform(rng) - 1.0;
+ q.y = 2.0*gsl_rng_uniform(rng) - 1.0;
+ q.z = 2.0*gsl_rng_uniform(rng) - 1.0;
+ q = normalise_quaternion(q);
+
+ return q;
+}
+
+
+/**
+ * quaternion_valid:
+ * @q: A %quaternion
+ *
+ * Checks if the given quaternion is normalised.
+ *
+ * This function performs a nasty floating point comparison of the form
+ * <code>(modulus > 0.999) && (modulus < 1.001)</code>, and so should not be
+ * relied upon to spot anything other than the most obvious input error.
+ *
+ * Returns: 1 if the quaternion is normalised, 0 if not.
+ **/
+int quaternion_valid(struct quaternion q)
+{
+ double qmod;
+
+ qmod = quaternion_modulus(q);
+
+ /* Modulus = 1 to within some tolerance?
+ * Nasty allowance for floating-point accuracy follows... */
+ if ( (qmod > 0.999) && (qmod < 1.001) ) return 1;
+
+ return 0;
+}
+
+
+/**
+ * quat_rot
+ * @q: A vector (in the form of a "struct rvec")
+ * @z: A %quaternion
+ *
+ * Rotates a vector according to a quaternion.
+ *
+ * Returns: A rotated version of @p.
+ **/
+struct rvec quat_rot(struct rvec q, struct quaternion z)
+{
+ struct rvec res;
+ double t01, t02, t03, t11, t12, t13, t22, t23, t33;
+
+ t01 = z.w*z.x;
+ t02 = z.w*z.y;
+ t03 = z.w*z.z;
+ t11 = z.x*z.x;
+ t12 = z.x*z.y;
+ t13 = z.x*z.z;
+ t22 = z.y*z.y;
+ t23 = z.y*z.z;
+ t33 = z.z*z.z;
+
+ res.u = (1.0 - 2.0 * (t22 + t33)) * q.u
+ + (2.0 * (t12 + t03)) * q.v
+ + (2.0 * (t13 - t02)) * q.w;
+
+ res.v = (2.0 * (t12 - t03)) * q.u
+ + (1.0 - 2.0 * (t11 + t33)) * q.v
+ + (2.0 * (t01 + t23)) * q.w;
+
+ res.w = (2.0 * (t02 + t13)) * q.u
+ + (2.0 * (t23 - t01)) * q.v
+ + (1.0 - 2.0 * (t11 + t22)) * q.w;
+
+ return res;
+}
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index 4955f875..a759ff15 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -61,44 +61,10 @@
#define THOMSON_LENGTH (2.81794e-15)
-/* ------------------------------ Quaternions ------------------------------- */
-
-/**
- * quaternion:
- *
- * <programlisting>
- * struct quaternion
- * {
- * double w
- * double x
- * double y
- * double z
- * };
- * </programlisting>
- *
- * A structure representing a quaternion.
- *
- **/
-struct quaternion;
-
-struct quaternion {
- double w;
- double x;
- double y;
- double z;
-};
-
#ifdef __cplusplus
extern "C" {
#endif
-extern struct quaternion normalise_quaternion(struct quaternion q);
-extern double quaternion_modulus(struct quaternion q);
-extern struct quaternion random_quaternion(gsl_rng *rng);
-extern int quaternion_valid(struct quaternion q);
-extern struct rvec quat_rot(struct rvec q, struct quaternion z);
-
-
/* --------------------------- Useful functions ----------------------------- */
extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v);
@@ -126,17 +92,6 @@ extern double flat_noise(gsl_rng *rng, double expected, double width);
extern double gaussian_noise(gsl_rng *rng, double expected, double stddev);
extern int poisson_noise(gsl_rng *rng, double expected);
-/* Keep these ones inline, to avoid function call overhead */
-static inline struct quaternion invalid_quaternion(void)
-{
- struct quaternion quat;
- quat.w = 0.0;
- quat.x = 0.0;
- quat.y = 0.0;
- quat.z = 0.0;
- return quat;
-}
-
static inline double distance(double x1, double y1, double x2, double y2)
{
return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
@@ -275,6 +230,47 @@ extern char *safe_basename(const char *in);
#define unlikely(x) (x)
#endif
+
+/* ------------------------------ Quaternions ------------------------------- */
+
+/**
+ * quaternion:
+ * @w: component
+ * @x: component
+ * @y: component
+ * @z: component
+ *
+ * A structure representing a quaternion.
+ *
+ **/
+struct quaternion;
+
+struct quaternion {
+ double w;
+ double x;
+ double y;
+ double z;
+};
+
+extern struct quaternion normalise_quaternion(struct quaternion q);
+extern double quaternion_modulus(struct quaternion q);
+extern struct quaternion random_quaternion(gsl_rng *rng);
+extern int quaternion_valid(struct quaternion q);
+extern struct rvec quat_rot(struct rvec q, struct quaternion z);
+
+/* Keep these ones inline, to avoid function call overhead */
+static inline struct quaternion invalid_quaternion(void)
+{
+ struct quaternion quat;
+ quat.w = 0.0;
+ quat.x = 0.0;
+ quat.y = 0.0;
+ quat.z = 0.0;
+ return quat;
+}
+
+
+
#ifdef __cplusplus
}
#endif
diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c
index 15826760..8a7daab7 100644
--- a/libcrystfel/src/xds.c
+++ b/libcrystfel/src/xds.c
@@ -3,12 +3,12 @@
*
* Invoke xds for crystal autoindexing
*
- * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2013 Cornelius Gati
*
* Authors:
- * 2010-2016 Thomas White <taw@physics.org>
+ * 2010-2017 Thomas White <taw@physics.org>
* 2013 Cornelius Gati <cornelius.gati@cfel.de>
*
* This file is part of CrystFEL.
@@ -498,7 +498,7 @@ static int write_inp(struct image *image, struct xds_private *xp)
}
-int run_xds(struct image *image, IndexingPrivate *priv)
+int run_xds(struct image *image, void *priv)
{
unsigned int opts;
int status;
@@ -619,8 +619,8 @@ int run_xds(struct image *image, IndexingPrivate *priv)
}
-IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell,
- struct detector *det, float *ltl)
+void *xds_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct detector *det, float *ltl)
{
struct xds_private *xp;
int need_cell = 0;
@@ -687,11 +687,11 @@ IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell,
xp->cell = cell;
xp->indm = *indm;
- return (IndexingPrivate *)xp;
+ return xp;
}
-void xds_cleanup(IndexingPrivate *pp)
+void xds_cleanup(void *pp)
{
struct xds_private *xp;
diff --git a/libcrystfel/src/xds.h b/libcrystfel/src/xds.h
index a0db2054..094d6d2f 100644
--- a/libcrystfel/src/xds.h
+++ b/libcrystfel/src/xds.h
@@ -3,13 +3,13 @@
*
* Invoke xds for crystal autoindexing
*
- * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
* Copyright © 2013 Cornelius Gati
*
* Authors:
- * 2010-2013 Thomas White <taw@physics.org>
- * 2013 Cornelius Gati <cornelius.gati@cfel.de>
+ * 2010-2013,2017 Thomas White <taw@physics.org>
+ * 2013 Cornelius Gati <cornelius.gati@cfel.de>
*
* This file is part of CrystFEL.
*
@@ -42,12 +42,12 @@
extern "C" {
#endif
-extern int run_xds(struct image *image, IndexingPrivate *ipriv);
+extern int run_xds(struct image *image, void *ipriv);
-extern IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell,
- struct detector *det, float *ltl);
+extern void *xds_prepare(IndexingMethod *indm, UnitCell *cell,
+ struct detector *det, float *ltl);
-extern void xds_cleanup(IndexingPrivate *pp);
+extern void xds_cleanup(void *pp);
#ifdef __cplusplus
}