aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/cell_explorer.c11
-rw-r--r--src/check_hkl.c9
-rw-r--r--src/cl-utils.c9
-rw-r--r--src/compare_hkl.c30
-rw-r--r--src/diffraction-gpu.c40
-rw-r--r--src/dw-hdfsee.c66
-rw-r--r--src/dw-hdfsee.h8
-rw-r--r--src/geoptimiser.c32
-rw-r--r--src/get_hkl.c4
-rw-r--r--src/im-sandbox.c4
-rw-r--r--src/indexamajig.c224
-rw-r--r--src/merge.c5
-rw-r--r--src/partialator.c6
-rw-r--r--src/pattern_sim.c2
-rw-r--r--src/process_hkl.c3
-rw-r--r--src/process_image.c77
-rw-r--r--src/process_image.h19
17 files changed, 379 insertions, 170 deletions
diff --git a/src/cell_explorer.c b/src/cell_explorer.c
index c8bbb069..f6c2c590 100644
--- a/src/cell_explorer.c
+++ b/src/cell_explorer.c
@@ -658,7 +658,9 @@ static void scan_cells(CellWindow *w)
if ( ignore ) continue;
- cell_get_parameters(cells[i], &a, &b, &c, &al, &be, &ga);
+ if ( cell_get_parameters(cells[i], &a, &b, &c, &al, &be, &ga) ) {
+ continue;
+ }
a *= 1e10; b *= 1e10; c *= 1e10;
al = rad2deg(al); be = rad2deg(be); ga = rad2deg(ga);
@@ -719,7 +721,10 @@ static void scan_minmax(CellWindow *w)
int j;
int found = 0;
- cell_get_parameters(w->cells[i], &a, &b, &c, &al, &be, &ga);
+ if ( cell_get_parameters(w->cells[i], &a, &b, &c, &al, &be, &ga) ) {
+ ERROR("Cell %i is bad\n", i);
+ continue;
+ }
a *= 1e10; b *= 1e10; c *= 1e10;
al = rad2deg(al); be = rad2deg(be); ga = rad2deg(ga);
@@ -1555,6 +1560,8 @@ int main(int argc, char *argv[])
return 1;
}
+ gsl_set_error_handler_off();
+
w.cells = NULL;
w.indms = NULL;
w.n_cells = 0;
diff --git a/src/check_hkl.c b/src/check_hkl.c
index def986c9..c7afb6cc 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -501,7 +501,8 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
signed int hs, ks, ls;
int bin;
- d = 2.0 * resolution(cell, h, k, l);
+ get_asymm(sym, h, k, l, &hs, &ks, &ls);
+ d = 2.0 * resolution(cell, hs, ks, ls);
if ( forbidden_reflection(cell, h, k, l) ) continue;
@@ -514,7 +515,6 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
}
if ( bin == -1 ) continue;
- get_asymm(sym, h, k, l, &hs, &ks, &ls);
if ( find_refl(counted, hs, ks, ls) != NULL ) continue;
add_refl(counted, hs, ks, ls);
@@ -530,7 +530,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
refl != NULL;
refl = next_refl(refl, iter) )
{
- signed int h, k, l;
+ signed int h, k, l, hs, ks, ls;
double d, val, esd;
int bin;
int j;
@@ -538,7 +538,8 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
get_indices(refl, &h, &k, &l);
if ( forbidden_reflection(cell, h, k, l) ) continue;
- d = resolution(cell, h, k, l) * 2.0;
+ get_asymm(sym, h, k, l, &hs, &ks, &ls);
+ d = resolution(cell, hs, ks, ls) * 2.0;
val = get_intensity(refl);
esd = get_esd_intensity(refl);
diff --git a/src/cl-utils.c b/src/cl-utils.c
index dc0a2713..a7e500cd 100644
--- a/src/cl-utils.c
+++ b/src/cl-utils.c
@@ -50,6 +50,15 @@ const char *clError(cl_int err)
case CL_SUCCESS :
return "no error";
+ case CL_DEVICE_NOT_AVAILABLE :
+ return "device not available";
+
+ case CL_DEVICE_NOT_FOUND :
+ return "device not found";
+
+ case CL_INVALID_DEVICE_TYPE :
+ return "invalid device type";
+
case CL_INVALID_PLATFORM :
return "invalid platform";
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 51e51dc9..81f12c94 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -286,13 +286,13 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2,
break;
case FOM_D1SIG :
- if ( fabs(i1-i2) < (sig1+sig2)/2.0 ) {
+ if ( fabs(i1-i2) < sqrt(sig1*sig1 + sig2*sig2) ) {
fctx->n_within[bin]++;
}
break;
case FOM_D2SIG :
- if ( fabs(i1-i2) < sig1+sig2 ) { /* = 2 * (sig1+sig2)/2 */
+ if ( fabs(i1-i2) < 2.0*sqrt(sig1*sig1 + sig2*sig2) ) {
fctx->n_within[bin]++;
}
break;
@@ -1125,7 +1125,7 @@ int main(int argc, char *argv[])
char *bfile = NULL;
char *sym_str = NULL;
SymOpList *sym;
- int ncom, nrej, nneg, nres, nbij, ncen;
+ int ncom, nrej, nmul, nneg, nres, nbij, ncen;
RefList *list1_acc;
RefList *list2_acc;
RefList *list1;
@@ -1149,6 +1149,7 @@ int main(int argc, char *argv[])
double min_I = +INFINITY;
double max_I = -INFINITY;
float highres, lowres;
+ int mul_cutoff = 0;
/* Long options */
const struct option longopts[] = {
@@ -1164,6 +1165,7 @@ int main(int argc, char *argv[])
{"shell-file", 1, NULL, 7},
{"highres", 1, NULL, 8},
{"lowres", 1, NULL, 9},
+ {"min-measurements", 1, NULL, 11},
{"ignore-negs", 0, &config_ignorenegs, 1},
{"zero-negs", 0, &config_zeronegs, 1},
{"intensity-shells", 0, &config_intshells, 1},
@@ -1260,6 +1262,13 @@ int main(int argc, char *argv[])
rmin_fix = 1.0 / (lowres/1e10);
break;
+ case 11 :
+ if ( sscanf(optarg, "%i", &mul_cutoff) != 1 ) {
+ ERROR("Invalid value for --min-measurements\n");
+ return 1;
+ }
+ break;
+
case '?' :
break;
@@ -1413,6 +1422,7 @@ int main(int argc, char *argv[])
/* Select reflections to be used */
ncom = 0;
nrej = 0;
+ nmul = 0;
nneg = 0;
nres = 0;
nbij = 0;
@@ -1426,6 +1436,7 @@ int main(int argc, char *argv[])
signed int h, k, l;
double val1, val2;
double esd1, esd2;
+ int mul1, mul2;
Reflection *refl2;
Reflection *refl1_acc;
Reflection *refl2_acc;
@@ -1441,6 +1452,9 @@ int main(int argc, char *argv[])
esd1 = get_esd_intensity(refl1);
esd2 = get_esd_intensity(refl2);
+ mul1 = get_redundancy(refl1);
+ mul2 = get_redundancy(refl2);
+
if ( (val1 < sigma_cutoff * esd1)
|| (val2 < sigma_cutoff * esd2) )
{
@@ -1453,6 +1467,11 @@ int main(int argc, char *argv[])
continue;
}
+ if ( (mul1 < mul_cutoff) || (mul2 < mul_cutoff) ) {
+ nmul++;
+ continue;
+ }
+
if ( config_zeronegs ) {
int d = 0;
if ( val1 < 0.0 ) {
@@ -1544,6 +1563,11 @@ int main(int argc, char *argv[])
" negative intensities which were set to zero.\n", nneg);
}
+ if ( nmul > 0 ) {
+ STATUS("%i reflection pairs rejected because either or both"
+ " versions had too few measurements.\n", nmul);
+ }
+
if ( nres > 0 ) {
STATUS("%i reflection pairs rejected because either or both"
" versions were outside the resolution range.\n", nres);
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 331170ae..22abcfd2 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -379,6 +379,7 @@ struct gpu_context *setup_gpu(int no_sfac,
float *flags_ptr;
size_t maxwgsize;
int i;
+ int have_ctx = 0;
char cflags[512] = "";
char *insert_stuff = NULL;
@@ -393,16 +394,37 @@ struct gpu_context *setup_gpu(int no_sfac,
ERROR("Couldn't find at least one platform!\n");
return NULL;
}
- prop[0] = CL_CONTEXT_PLATFORM;
- prop[1] = (cl_context_properties)platforms[0];
- prop[2] = 0;
- gctx = malloc(sizeof(*gctx));
- gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU,
- NULL, NULL, &err);
- if ( err != CL_SUCCESS ) {
- ERROR("Couldn't create OpenCL context: %i\n", err);
- free(gctx);
+ /* Find a GPU platform in the list */
+ for ( i=0; i<nplat; i++ ) {
+
+ prop[0] = CL_CONTEXT_PLATFORM;
+ prop[1] = (cl_context_properties)platforms[i];
+ prop[2] = 0;
+
+ gctx = malloc(sizeof(*gctx));
+ gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU,
+ NULL, NULL, &err);
+
+ if ( err != CL_SUCCESS ) {
+ if ( err == CL_DEVICE_NOT_FOUND ) {
+ /* No GPU device in this platform */
+ continue; /* Try next platform */
+ } else {
+ ERROR("Couldn't create OpenCL context: %s (%i)\n",
+ clError(err), err);
+ free(gctx);
+ return NULL;
+ }
+ } else {
+ STATUS("Using OpenCL platform %i (%i total)\n", i, nplat);
+ have_ctx = 1;
+ break;
+ }
+ }
+
+ if ( !have_ctx ) {
+ ERROR("Couldn't find a GPU device in any platform.\n");
return NULL;
}
diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c
index 87c23f4e..0336f25e 100644
--- a/src/dw-hdfsee.c
+++ b/src/dw-hdfsee.c
@@ -74,8 +74,8 @@ static void displaywindow_error(DisplayWindow *dw, const char *message)
static gint displaywindow_closed(GtkWidget *window, DisplayWindow *dw)
{
- if ( dw->hdfile != NULL ) {
- hdfile_close(dw->hdfile);
+ if ( dw->imagefile != NULL ) {
+ imagefile_close(dw->imagefile);
}
if ( dw->surf != NULL ) cairo_surface_destroy(dw->surf);
@@ -725,7 +725,7 @@ static gint displaywindow_set_binning(GtkWidget *widget, DisplayWindow *dw)
return 0;
}
- if ( dw->hdfile == NULL ) {
+ if ( dw->imagefile == NULL ) {
return 0;
}
@@ -852,8 +852,8 @@ static gint displaywindow_set_boostint(GtkWidget *widget, DisplayWindow *dw)
return 0;
}
- if ( dw->hdfile == NULL ) {
- return 0;
+ if ( dw->imagefile == NULL ) {
+ return 0;
}
bd = malloc(sizeof(BoostIntDialog));
@@ -927,8 +927,8 @@ static gint displaywindow_newevent(DisplayWindow *dw, int new_event)
float **old_dp = dw->image->dp;
int **old_bad = dw->image->bad;
- fail = hdf5_read2(dw->hdfile, dw->image,
- dw->ev_list->events[new_event], 0);
+ fail = imagefile_read(dw->imagefile, dw->image,
+ dw->ev_list->events[new_event]);
if ( fail ) {
ERROR("Couldn't load image");
dw->image->dp = old_dp;
@@ -1046,8 +1046,8 @@ static gint displaywindow_set_newevent(GtkWidget *widget, DisplayWindow *dw)
return 0;
}
- if ( dw->hdfile == NULL ) {
- return 0;
+ if ( dw->imagefile == NULL ) {
+ return 0;
}
ed = malloc(sizeof(EventDialog));
@@ -1162,7 +1162,7 @@ static gint displaywindow_set_ringradius(GtkWidget *widget, DisplayWindow *dw)
return 0;
}
- if ( dw->hdfile == NULL ) {
+ if ( dw->imagefile == NULL ) {
return 0;
}
@@ -1796,7 +1796,7 @@ static gint displaywindow_show_numbers(GtkWidget *widget, DisplayWindow *dw)
return 0;
}
- if ( dw->hdfile == NULL ) {
+ if ( dw->imagefile == NULL ) {
return 0;
}
@@ -2306,12 +2306,14 @@ struct newhdf {
char name[1024];
};
+/* New HDF5 element selected from menu */
static gint displaywindow_newhdf(GtkMenuItem *item, struct newhdf *nh)
{
gboolean a;
int fail;
if ( nh->dw->not_ready_yet ) return 0;
+ assert(imagefile_get_type(nh->dw->imagefile) == IMAGEFILE_HDF5);
a = gtk_check_menu_item_get_active(GTK_CHECK_MENU_ITEM(nh->widget));
if ( !a ) return 0;
@@ -2320,7 +2322,8 @@ static gint displaywindow_newhdf(GtkMenuItem *item, struct newhdf *nh)
* one */
free_detector_geometry(nh->dw->image->det);
nh->dw->image->det = NULL;
- fail = hdf5_read(nh->dw->hdfile, nh->dw->image, nh->name, 0);
+ fail = hdf5_read(imagefile_get_hdfile(nh->dw->imagefile),
+ nh->dw->image, nh->name, 0);
if ( fail ) {
ERROR("Couldn't load image");
return 1;
@@ -2456,7 +2459,12 @@ static int displaywindow_update_menus(DisplayWindow *dw, const char *selectme)
GtkWidget *ms;
GtkWidget *w;
- ms = displaywindow_createhdfmenus(dw->hdfile, dw, selectme);
+ if ( imagefile_get_type(dw->imagefile) != IMAGEFILE_HDF5 ) {
+ return 0;
+ }
+
+ ms = displaywindow_createhdfmenus(imagefile_get_hdfile(dw->imagefile),
+ dw, selectme);
if ( ms == NULL ) return 1;
@@ -2828,8 +2836,8 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename,
dw->image->lambda = 0.0;
dw->image->filename = filename;
- dw->hdfile = hdfile_open(filename);
- if ( dw->hdfile == NULL ) {
+ dw->imagefile = imagefile_open(filename);
+ if ( dw->imagefile == NULL ) {
ERROR("Couldn't open file: %s\n", filename);
free(dw->geom_filename);
free(dw);
@@ -2842,11 +2850,19 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename,
}
if ( dw->image->det != NULL && ( dw->image->det->path_dim != 0 ||
- dw->image->det->dim_dim != 0 )) {
+ dw->image->det->dim_dim != 0 ))
+ {
+ struct hdfile *hdfile;
+
+ if ( imagefile_get_type(dw->imagefile) != IMAGEFILE_HDF5 ) {
+ ERROR("Multi-event geometry, but not HDF5 file!\n");
+ return NULL;
+ }
+ hdfile = imagefile_get_hdfile(dw->imagefile);
dw->multi_event = 1;
- dw->ev_list = fill_event_list(dw->hdfile, dw->image->det);
+ dw->ev_list = fill_event_list(hdfile, dw->image->det);
if ( dw->ev_list == NULL ) {
ERROR("Error while parsing file structure\n");
@@ -2883,18 +2899,26 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename,
dw->curr_event = 0;
ev = dw->ev_list->events[dw->curr_event];
}
- check = hdf5_read2(dw->hdfile, dw->image, ev, 0);
+ check = imagefile_read(dw->imagefile, dw->image, ev);
} else {
- check = hdf5_read2(dw->hdfile, dw->image, NULL, 0);
+ check = imagefile_read(dw->imagefile, dw->image, NULL);
}
} else {
- check = hdf5_read(dw->hdfile, dw->image, element, 0);
+ if ( element != NULL ) {
+ if ( imagefile_get_type(dw->imagefile) != IMAGEFILE_HDF5 ) {
+ ERROR("-e/--image requiest an HDF5 file\n");
+ return NULL;
+ }
+ hdfile_set_image(imagefile_get_hdfile(dw->imagefile),
+ element);
+ }
+ check = imagefile_read_simple(dw->imagefile, dw->image);
dw->simple = 1;
}
if ( check ) {
ERROR("Couldn't load file\n");
- hdfile_close(dw->hdfile);
+ imagefile_close(dw->imagefile);
free(dw->geom_filename);
return NULL;
}
diff --git a/src/dw-hdfsee.h b/src/dw-hdfsee.h
index 85c82ac7..fdbea231 100644
--- a/src/dw-hdfsee.h
+++ b/src/dw-hdfsee.h
@@ -3,12 +3,12 @@
*
* Quick yet non-crappy HDF viewer
*
- * 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:
- * 2009-2014 Thomas White <taw@physics.org>
+ * 2009-2017 Thomas White <taw@physics.org>
* 2014 Valerio Mariani
* 2014 Takanori Nakane <nakane.t@gmail.com>
* 2012 Richard Kirian
@@ -101,8 +101,8 @@ typedef struct {
int simple;
- struct hdfile *hdfile;
- struct image *image;
+ struct imagefile *imagefile;
+ struct image *image;
char *geom_filename;
char *rg_coll_name;
diff --git a/src/geoptimiser.c b/src/geoptimiser.c
index e0337358..da9b16c4 100644
--- a/src/geoptimiser.c
+++ b/src/geoptimiser.c
@@ -631,10 +631,18 @@ static int add_distance_to_list(struct gpanel *gp,
double *det_shift)
{
int pix_index;
+ int ifs, iss;
double rfs, rss;
double crx, cry;
- pix_index = ((int)rint(imfe->fs) + gp->p->w*(int)rint(imfe->ss));
+ ifs = imfe->fs;
+ iss = imfe->ss; /* Explicit rounding towards zero (truncation) */
+ pix_index = ifs + gp->p->w*iss;
+
+ if ( (ifs >= gp->p->w) || (iss >= gp->p->h) ) {
+ ERROR("Peak is outside panel!\n");
+ return 1;
+ }
if ( gp->num_pix_displ[pix_index] > 0 ) {
@@ -782,7 +790,14 @@ static int compute_pixel_displacements(struct image *images, int n_images,
r = add_distance_to_list(gp, imfe, refl, fx, fy,
det_shift);
- if ( r ) return r;
+ if ( r ) {
+ ERROR("Error processing peak %f,%f "
+ "(panel %s), image %s %s\n",
+ imfe->fs, imfe->ss, gp->p->name,
+ images[cp].filename,
+ get_event_string(images[cp].event));
+ return r;
+ }
}
}
@@ -1070,12 +1085,7 @@ static void correct_rotation_and_stretch(struct rg_collection *connected,
/* Calculate corner adjustment
* NB All panels follow the first one */
- if ( ip == 0 ) {
-
- new_cnx = p->cnx * cs;
- new_cny = p->cny * cs;
-
- } else {
+ if ( ip > 0 && connected->rigid_groups[di]->n_panels == 2 && !gparams->no_cspad ) {
struct panel *p0;
double delta_x, delta_y, delta;
@@ -1090,6 +1100,10 @@ static void correct_rotation_and_stretch(struct rg_collection *connected,
new_cnx = p0->cnx + delta*p0->fsx;
new_cny = p0->cny + delta*p0->fsy;
+ } else {
+
+ new_cnx = p->cnx * cs;
+ new_cny = p->cny * cs;
}
/* The average displacements now need to be updated */
@@ -2488,8 +2502,10 @@ int main(int argc, char *argv[])
}
#ifdef HAVE_SAVE_TO_PNG
+#if !GLIB_CHECK_VERSION(2,36,0)
g_type_init();
#endif
+#endif
ret_val = optimize_geometry(gparams, det, quadrants, connected);
diff --git a/src/get_hkl.c b/src/get_hkl.c
index d6efe747..77f34da2 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -98,7 +98,9 @@ static void show_help(const char *s)
static void copy_notes(RefList *out, RefList *in)
{
- reflist_add_notes(out, reflist_get_notes(in));
+ if ( reflist_get_notes(in) != NULL ) {
+ reflist_add_notes(out, reflist_get_notes(in));
+ }
}
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index f5493453..60400ad8 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -372,10 +372,10 @@ static void run_work(const struct index_args *iargs, Stream *st,
}
time_accounts_set(taccs, TACC_FINALCLEANUP);
- cleanup_indexing(iargs->indm, iargs->ipriv);
+ cleanup_indexing(iargs->ipriv);
free_detector_geometry(iargs->det);
free(iargs->hdf5_peak_path);
- free_copy_hdf5_field_list(iargs->copyme);
+ free_imagefile_field_list(iargs->copyme);
cell_free(iargs->cell);
if ( iargs->profile ) time_accounts_print(taccs);
time_accounts_free(taccs);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index a1997b9e..ce87fde6 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -3,16 +3,17 @@
*
* Index patterns, output hkl+intensity etc.
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
* Copyright © 2012 Lorenzo Galli
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2017 Thomas White <taw@physics.org>
* 2011 Richard Kirian
* 2012 Lorenzo Galli
* 2012 Chunhong Yoon
+ * 2017 Valerio Mariani <valerio.mariani@desy.de>
*
* This file is part of CrystFEL.
*
@@ -88,6 +89,7 @@ static void show_help(const char *s)
" --peaks=<method> Use 'method' for finding peaks. Choose from:\n"
" zaef : Use Zaefferer (2000) gradient detection.\n"
" This is the default method.\n"
+" peakfinder8: Use Peakfinder8 algorithm.\n"
" hdf5 : Get from a table in HDF5 file.\n"
" cxi : Get from CXI format HDF5 file.\n"
" --hdf5-peaks=<p> Find peaks table in HDF5 file here.\n"
@@ -95,32 +97,52 @@ static void show_help(const char *s)
" --integration=<meth> Perform final pattern integration using <meth>.\n"
"\n\n"
"For more control over the process, you might need:\n\n"
-" --tolerance=<tol> Set the tolerances for cell comparison.\n"
-" Default: 5,5,5,1.5.\n"
-" --filter-noise Apply an aggressive noise filter which sets all\n"
-" pixels in each 3x3 region to zero if any of them\n"
-" have negative values. Intensity measurement will\n"
-" be performed on the image as it was before this.\n"
-" --median-filter=<n> Apply a median filter to the image data. Intensity\n"
-" measurement will be performed on the image as it\n"
-" was before this. The side length of the median\n"
-" filter box will be 2<n>+1. Default: 0 (no filter).\n"
-" --no-sat-corr Don't correct values of saturated peaks using a\n"
-" table included in the HDF5 file.\n"
-" --threshold=<n> Only accept peaks above <n> ADU. Default: 800.\n"
-" --min-gradient=<n> Minimum squared gradient for Zaefferer peak search.\n"
-" Default: 100,000.\n"
-" --min-snr=<n> Minimum signal-to-noise ratio for peaks.\n"
-" Default: 5.\n"
-" --check-hdf5-snr Check SNR for peaks from --peaks=hdf5.\n"
-" --peak-radius=<r> Integration radii for peak search.\n"
-" --int-radius=<r> Set the integration radii. Default: 4,5,7.\n"
-" --push-res=<n> Integrate higher than apparent resolution cutoff.\n"
-" --highres=<n> Absolute resolution cutoff in Angstroms.\n"
-" --fix-profile-radius Fix the reciprocal space profile radius for spot\n"
-" prediction (default: automatically determine.\n"
-" --fix-bandwidth Set the bandwidth for spot prediction.\n"
-" --fix-divergence Set the divergence (full angle) for spot prediction.\n"
+" --tolerance=<tol> Set the tolerances for cell comparison.\n"
+" Default: 5,5,5,1.5.\n"
+" --filter-noise Apply an aggressive noise filter which sets all\n"
+" pixels in each 3x3 region to zero if any of them\n"
+" have negative values. Intensity measurement will\n"
+" be performed on the image as it was before this.\n"
+" --median-filter=<n> Apply a median filter to the image data. Intensity\n"
+" measurement will be performed on the image as it\n"
+" was before this. The side length of the median\n"
+" filter box will be 2<n>+1.\n"
+" Default: 0 (no filter).\n"
+" --no-sat-corr Don't correct values of saturated peaks using a\n"
+" table included in the HDF5 file.\n"
+" --threshold=<n> Only accept peaks above <n> ADU in both the\n"
+" Zaefferer and Peakfinder8 algorithms.\n"
+" Default: 800.\n"
+" --min-gradient=<n> Minimum squared gradient for Zaefferer peak\n"
+" search. Default: 100,000.\n"
+" --min-snr=<n> Minimum signal-to-noise ratio for peaks, with both\n"
+" Zaefferer and Peakfinder8 algorithms.\n"
+" Default: 5.\n"
+" --min-pix-count=<n> Only accept peaks if they include more than <n>\n"
+" pixels, in the Peakfinder8 algorithm.\n"
+" Default: 2.\n"
+" --max-pix-count=<n> Only accept peaks if they include less than <n>\n"
+" pixels, in the Peakfinder8 algorithm.\n"
+" Default: 200.\n"
+" --min-peaks=<n> Minimum number of peaks for indexing.\n"
+" --local-bg-radius=<n> Radius (in pixels) to use for the estimation of\n"
+" local background in the Peakfinder8 algorithm.\n"
+" Default: 3.\n"
+" --min-res=<n> Only accept peaks if they lay at more than <n>\n"
+" pixels from the center of the detector, in the\n"
+" peakfinder8 algorithm. Default: 0.\n"
+" --max-res=<n> Only accept peaks if they lay at less than <n>\n"
+" pixels from the center of the detector, in the\n"
+" peakfinder8 algorithm. Default: 1200.\n"
+" --check-hdf5-snr Check SNR for peaks from --peaks=hdf5.\n"
+" --peak-radius=<r> Integration radii for peak search.\n"
+" --int-radius=<r> Set the integration radii. Default: 4,5,7.\n"
+" --push-res=<n> Integrate higher than apparent resolution cutoff.\n"
+" --highres=<n> Absolute resolution cutoff in Angstroms.\n"
+" --fix-profile-radius Fix the reciprocal space profile radius for spot\n"
+" prediction (default: automatically determine.\n"
+" --fix-bandwidth Set the bandwidth for spot prediction.\n"
+" --fix-divergence Set the divergence (full angle) for spot prediction.\n"
"\n"
"\nFor time-resolved stuff, you might want to use:\n\n"
" --copy-hdf5-field <f> Copy the value of field <f> into the stream. You\n"
@@ -131,16 +153,19 @@ static void show_help(const char *s)
" --temp-dir=<path> Put the temporary folder under <path>.\n"
"\n"
"\nOptions you probably won't need:\n\n"
-" --no-check-prefix Don't attempt to correct the --prefix.\n"
-" --no-use-saturated During the initial peak search, reject\n"
-" peaks which contain pixels above max_adu.\n"
-" --no-revalidate Don't re-integrate and check HDF5 peaks for\n"
-" validity.\n"
-" --no-peaks-in-stream Do not record peak search results in the stream.\n"
-" --no-refls-in-stream Do not record integrated reflections in the stream.\n"
-" --int-diag=<cond> Show debugging information about reflections.\n"
-" --no-refine Skip the prediction refinement step.\n"
-" --profile Show timing data for performance monitoring.\n"
+" --no-check-prefix Don't attempt to correct the --prefix.\n"
+" --no-use-saturated During the initial peak search, reject\n"
+" peaks which contain pixels above max_adu.\n"
+" --no-revalidate Don't re-integrate and check HDF5 peaks for\n"
+" validity.\n"
+" --no-peaks-in-stream Do not record peak search results in the stream.\n"
+" --no-refls-in-stream Do not record integrated reflections in the stream.\n"
+" --no-non-hits-in-stream Do not include non-hit frames in the stream.\n"
+" (see --min-peaks)\n"
+" --int-diag=<cond> Show debugging information about reflections.\n"
+" --no-refine Skip the prediction refinement step.\n"
+" --profile Show timing data for performance monitoring.\n"
+" --no-half-pixel-shift Don't offset the HDF5 peak locations by 0.5 px.\n"
"\nLow-level options for the felix indexer:\n\n"
" --felix-options Change the default arguments passed to the indexer.\n"
" Given as a list of comma separated list of \n"
@@ -150,9 +175,9 @@ static void show_help(const char *s)
}
-static void add_geom_beam_stuff_to_copy_hdf5(struct copy_hdf5_field *copyme,
- struct detector *det,
- struct beam_params *beam)
+static void add_geom_beam_stuff_to_field_list(struct imagefile_field_list *copyme,
+ struct detector *det,
+ struct beam_params *beam)
{
int i;
@@ -161,12 +186,12 @@ static void add_geom_beam_stuff_to_copy_hdf5(struct copy_hdf5_field *copyme,
struct panel *p = &det->panels[i];
if ( p->clen_from != NULL ) {
- add_copy_hdf5_field(copyme, p->clen_from);
+ add_imagefile_field(copyme, p->clen_from);
}
}
if ( beam->photon_energy_from != NULL ) {
- add_copy_hdf5_field(copyme, beam->photon_energy_from);
+ add_imagefile_field(copyme, beam->photon_energy_from);
}
}
@@ -181,8 +206,6 @@ int main(int argc, char *argv[])
int config_checkprefix = 1;
int config_basename = 0;
int integrate_saturated = 0;
- IndexingMethod *indm;
- IndexingPrivate **ipriv;
char *indm_str = NULL;
char *cellfile = NULL;
char *prefix = NULL;
@@ -214,11 +237,17 @@ int main(int argc, char *argv[])
iargs.threshold = 800.0;
iargs.min_gradient = 100000.0;
iargs.min_snr = 5.0;
+ iargs.min_pix_count = 2;
+ iargs.max_pix_count = 200;
+ iargs.min_res = 0;
+ iargs.max_res = 1200;
+ iargs.local_bg_radius = 3;
iargs.check_hdf5_snr = 0;
iargs.det = NULL;
iargs.peaks = PEAK_ZAEF;
iargs.beam = &beam;
iargs.hdf5_peak_path = NULL;
+ iargs.half_pixel_shift = 1;
iargs.copyme = NULL;
iargs.pk_inn = -1.0;
iargs.pk_mid = -1.0;
@@ -230,13 +259,14 @@ int main(int argc, char *argv[])
iargs.no_revalidate = 0;
iargs.stream_peaks = 1;
iargs.stream_refls = 1;
+ iargs.stream_nonhits = 1;
iargs.int_diag = INTDIAG_NONE;
- iargs.copyme = new_copy_hdf5_field_list();
+ iargs.copyme = new_imagefile_field_list();
+ iargs.min_peaks = 0;
if ( iargs.copyme == NULL ) {
ERROR("Couldn't allocate HDF5 field list.\n");
return 1;
}
- iargs.indm = NULL; /* No default */
iargs.ipriv = NULL; /* No default */
iargs.int_meth = integration_method("rings-nocen-nosat-nograd", NULL);
iargs.push_res = 0.0;
@@ -268,12 +298,14 @@ int main(int argc, char *argv[])
{"basename", 0, &config_basename, 1},
{"no-peaks-in-stream", 0, &iargs.stream_peaks, 0},
{"no-refls-in-stream", 0, &iargs.stream_refls, 0},
+ {"no-non-hits-in-stream", 0, &iargs.stream_nonhits, 0},
{"integrate-saturated",0, &integrate_saturated, 1},
{"no-use-saturated", 0, &iargs.use_saturated, 0},
{"no-revalidate", 0, &iargs.no_revalidate, 1},
{"check-hdf5-snr", 0, &iargs.check_hdf5_snr, 1},
{"no-refine", 0, &no_refine, 1},
{"profile", 0, &iargs.profile, 1},
+ {"no-half-pixel-shift",0, &iargs.half_pixel_shift, 0},
/* Long-only options which don't actually do anything */
{"no-sat-corr", 0, &iargs.satcorr, 0},
@@ -306,6 +338,12 @@ int main(int argc, char *argv[])
{"fix-bandwidth", 1, NULL, 23},
{"fix-divergence", 1, NULL, 24},
{"felix-options", 1, NULL, 25},
+ {"min-pix-count", 1, NULL, 26},
+ {"max-pix-count", 1, NULL, 27},
+ {"local-bg-radius", 1, NULL, 28},
+ {"min-res", 1, NULL, 29},
+ {"max-res", 1, NULL, 30},
+ {"min-peaks", 1, NULL, 31},
{0, 0, NULL, 0}
};
@@ -400,7 +438,7 @@ int main(int argc, char *argv[])
break;
case 10 :
- add_copy_hdf5_field(iargs.copyme, optarg);
+ add_imagefile_field(iargs.copyme, optarg);
break;
case 11 :
@@ -489,6 +527,30 @@ int main(int argc, char *argv[])
}
break;
+ case 26:
+ iargs.min_pix_count = atoi(optarg);
+ break;
+
+ case 27:
+ iargs.max_pix_count = atoi(optarg);
+ break;
+
+ case 28:
+ iargs.local_bg_radius = atoi(optarg);
+ break;
+
+ case 29:
+ iargs.min_res = atoi(optarg);
+ break;
+
+ case 30:
+ iargs.max_res = atoi(optarg);
+ break;
+
+ case 31:
+ iargs.min_peaks = atoi(optarg);
+ break;
+
case 0 :
break;
@@ -541,6 +603,8 @@ int main(int argc, char *argv[])
}
if ( strcmp(speaks, "zaef") == 0 ) {
iargs.peaks = PEAK_ZAEF;
+ } else if ( strcmp(speaks, "peakfinder8") == 0 ) {
+ iargs.peaks = PEAK_PEAKFINDER8;
} else if ( strcmp(speaks, "hdf5") == 0 ) {
iargs.peaks = PEAK_HDF5;
} else if ( strcmp(speaks, "cxi") == 0 ) {
@@ -574,7 +638,7 @@ int main(int argc, char *argv[])
geom_filename);
return 1;
}
- add_geom_beam_stuff_to_copy_hdf5(iargs.copyme, iargs.det, iargs.beam);
+ add_geom_beam_stuff_to_field_list(iargs.copyme, iargs.det, iargs.beam);
/* If no peak path from geometry file, use these (but see later) */
if ( iargs.hdf5_peak_path == NULL ) {
@@ -591,35 +655,6 @@ int main(int argc, char *argv[])
iargs.hdf5_peak_path = command_line_peak_path;
}
- /* Parse indexing methods */
- if ( indm_str == NULL ) {
-
- STATUS("You didn't specify an indexing method, so I won't try "
- " to index anything.\n"
- "If that isn't what you wanted, re-run with"
- " --indexing=<methods>.\n");
- indm = NULL;
-
- } else {
-
- int i = 0;
-
- indm = build_indexer_list(indm_str);
- if ( indm == NULL ) {
- ERROR("Invalid indexer list '%s'\n", indm_str);
- return 1;
- }
- free(indm_str);
-
- /* If --no-refine, unset the refinement flag on all methods */
- if ( no_refine ) {
- while ( indm[i] != INDEXING_NONE ) {
- indm[i] &= ~INDEXING_REFINE;
- i++;
- }
- }
- }
-
/* Parse integration method */
if ( int_str != NULL ) {
@@ -755,27 +790,34 @@ int main(int argc, char *argv[])
}
free(outfile);
- /* Prepare the indexer */
- if ( indm != NULL ) {
- ipriv = prepare_indexing(indm, iargs.cell, iargs.det,
- iargs.tols, iargs.felix_options);
- if ( ipriv == NULL ) {
- ERROR("Failed to prepare indexing.\n");
+ /* Prepare the indexing system */
+ if ( indm_str == NULL ) {
+
+ STATUS("You didn't specify an indexing method, so I won't try "
+ " to index anything.\n"
+ "If that isn't what you wanted, re-run with"
+ " --indexing=<methods>.\n");
+ iargs.ipriv = NULL;
+
+ } else {
+
+ iargs.ipriv = setup_indexing(indm_str, iargs.cell, iargs.det,
+ iargs.tols, no_refine,
+ iargs.felix_options);
+ if ( iargs.ipriv == NULL ) {
+ ERROR("Failed to set up indexing system\n");
return 1;
}
- } else {
- ipriv = NULL;
+ free(indm_str);
+
}
gsl_set_error_handler_off();
- iargs.indm = indm;
- iargs.ipriv = ipriv;
-
create_sandbox(&iargs, n_proc, prefix, config_basename, fh,
st, tempdir);
- free_copy_hdf5_field_list(iargs.copyme);
+ free_imagefile_field_list(iargs.copyme);
cell_free(iargs.cell);
free(iargs.beam->photon_energy_from);
free(prefix);
@@ -783,7 +825,7 @@ int main(int argc, char *argv[])
free_detector_geometry(iargs.det);
free(iargs.hdf5_peak_path);
close_stream(st);
- cleanup_indexing(indm, ipriv);
+ cleanup_indexing(iargs.ipriv);
return 0;
}
diff --git a/src/merge.c b/src/merge.c
index 8d1fae0f..9734c469 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -237,10 +237,7 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
Reflection *refl;
RefListIterator *iter;
- if ( n == 0 ) {
- ERROR("No crystals!\n");
- return NULL;
- }
+ if ( n == 0 ) return NULL;
full = reflist_new();
diff --git a/src/partialator.c b/src/partialator.c
index 09feebb4..569145e8 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -191,6 +191,12 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
snprintf(tmp, 1024, "%s1", outfile);
split = merge_intensities(crystals1, n_crystals1, nthreads,
pmodel, min_measurements, push_res, 1);
+
+ if ( split == NULL ) {
+ ERROR("Not enough crystals for two way split!\n");
+ return;
+ }
+
STATUS("Writing two-way split to %s ", tmp);
write_reflist_2(tmp, split, sym);
reflist_free(split);
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index e1e04789..e9fc9916 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -745,7 +745,7 @@ int main(int argc, char *argv[])
"each of a, b and c\n");
}
if ( intfile == NULL ) {
- STATUS(" Full intensities: all equal");
+ STATUS(" Full intensities: all equal\n");
} else {
STATUS(" Full intensities: from %s\n", intfile);
}
diff --git a/src/process_hkl.c b/src/process_hkl.c
index d2cf640b..301bc6e4 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -287,7 +287,8 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
if ( isnan(scale) ) return 1;
if ( scale <= 0.0 ) return 1;
if ( stat != NULL ) {
- fprintf(stat, "%s %f %f\n", image->filename, scale, cc);
+ fprintf(stat, "%s %s %f %f\n", image->filename,
+ get_event_string(image->event), scale, cc);
}
} else {
diff --git a/src/process_image.c b/src/process_image.c
index bcaee543..498b3398 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -8,7 +8,7 @@
*
* Authors:
* 2010-2017 Thomas White <taw@physics.org>
- * 2014 Valerio Mariani
+ * 2014-2017 Valerio Mariani <valerio.mariani@desy.de>
*
* This file is part of CrystFEL.
*
@@ -102,7 +102,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
int serial, struct sb_shm *sb_shared, TimeAccounts *taccs)
{
int check;
- struct hdfile *hdfile;
+ struct imagefile *imfile;
struct image image;
int i;
int r;
@@ -125,15 +125,15 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
time_accounts_set(taccs, TACC_HDF5OPEN);
sb_shared->pings[cookie]++;
- hdfile = hdfile_open(image.filename);
- if ( hdfile == NULL ) {
+ imfile = imagefile_open(image.filename);
+ if ( imfile == NULL ) {
ERROR("Couldn't open file: %s\n", image.filename);
return;
}
time_accounts_set(taccs, TACC_HDF5READ);
sb_shared->pings[cookie]++;
- check = hdf5_read2(hdfile, &image, image.event, 0);
+ check = imagefile_read(imfile, &image, image.event);
if ( check ) {
return;
}
@@ -159,8 +159,14 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
sb_shared->pings[cookie]++;
switch ( iargs->peaks ) {
+ struct hdfile *hdfile;
+
case PEAK_HDF5:
- if ( get_peaks(&image, hdfile, iargs->hdf5_peak_path) ) {
+ hdfile = imagefile_get_hdfile(imfile);
+ if ( (hdfile == NULL)
+ || (get_peaks_2(&image, hdfile, iargs->hdf5_peak_path,
+ iargs->half_pixel_shift)) )
+ {
ERROR("Failed to get peaks from HDF5 file.\n");
}
if ( !iargs->no_revalidate ) {
@@ -172,8 +178,12 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
break;
case PEAK_CXI:
- if ( get_peaks_cxi(&image, hdfile, iargs->hdf5_peak_path,
- pargs->filename_p_e) ) {
+ hdfile = imagefile_get_hdfile(imfile);
+ if ( (hdfile == NULL)
+ || (get_peaks_cxi_2(&image, hdfile, iargs->hdf5_peak_path,
+ pargs->filename_p_e,
+ iargs->half_pixel_shift)) )
+ {
ERROR("Failed to get peaks from CXI file.\n");
}
if ( !iargs->no_revalidate ) {
@@ -191,6 +201,28 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
iargs->use_saturated);
break;
+ case PEAK_PEAKFINDER8:
+ if ( search_peaks_peakfinder8(&image, 2048,
+ iargs->threshold,
+ iargs->min_snr,
+ iargs->min_pix_count,
+ iargs->max_pix_count,
+ iargs->local_bg_radius,
+ iargs->min_res,
+ iargs->max_res,
+ iargs->use_saturated) ) {
+ if ( image.event != NULL ) {
+ ERROR("Failed to find peaks in image %s"
+ "(event %s).\n", image.filename,
+ get_event_string(image.event));
+ } else {
+ ERROR("Failed to find peaks in image %s.",
+ image.filename);
+ }
+
+ }
+ break;
+
}
restore_image_data(image.dp, image.det, prefilter);
@@ -201,7 +233,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
if ( r ) {
ERROR("Failed to chdir to temporary folder: %s\n",
strerror(errno));
- hdfile_close(hdfile);
+ imagefile_close(imfile);
return;
}
@@ -217,15 +249,30 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
image.bw = 0.00000001;
}
+ if ( image_feature_count(image.features) < iargs->min_peaks ) {
+ r = chdir(rn);
+ if ( r ) {
+ ERROR("Failed to chdir: %s\n", strerror(errno));
+ imagefile_close(imfile);
+ return;
+ }
+ free(rn);
+
+ if ( iargs->stream_nonhits ) {
+ goto streamwrite;
+ } else {
+ goto out;
+ }
+ }
+
/* Index the pattern */
time_accounts_set(taccs, TACC_INDEXING);
- index_pattern_2(&image, iargs->indm, iargs->ipriv,
- &sb_shared->pings[cookie]);
+ index_pattern_2(&image, iargs->ipriv, &sb_shared->pings[cookie]);
r = chdir(rn);
if ( r ) {
ERROR("Failed to chdir: %s\n", strerror(errno));
- hdfile_close(hdfile);
+ imagefile_close(imfile);
return;
}
free(rn);
@@ -261,9 +308,10 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
iargs->int_diag_k, iargs->int_diag_l,
&sb_shared->term_lock);
+streamwrite:
time_accounts_set(taccs, TACC_WRITESTREAM);
sb_shared->pings[cookie]++;
- ret = write_chunk(st, &image, hdfile,
+ ret = write_chunk(st, &image, imfile,
iargs->stream_peaks, iargs->stream_refls,
pargs->filename_p_e->ev);
if ( ret != 0 ) {
@@ -280,6 +328,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
get_event_string(image.event));
}
+out:
/* Count crystals which are still good */
time_accounts_set(taccs, TACC_TOTALS);
sb_shared->pings[cookie]++;
@@ -313,5 +362,5 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
image_feature_list_free(image.features);
free_detector_geometry(image.det);
- hdfile_close(hdfile);
+ imagefile_close(imfile);
}
diff --git a/src/process_image.h b/src/process_image.h
index d41c23f5..ec51c188 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -3,12 +3,12 @@
*
* The processing pipeline for one image
*
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
* 2010-2016 Thomas White <taw@physics.org>
- * 2014 Valerio Mariani
+ * 2014-2017 Valerio Mariani <valerio.mariani@desy.de>
*
* This file is part of CrystFEL.
*
@@ -42,6 +42,7 @@ struct index_args;
enum {
+ PEAK_PEAKFINDER8,
PEAK_ZAEF,
PEAK_HDF5,
PEAK_CXI,
@@ -61,24 +62,32 @@ struct index_args
float min_snr;
int check_hdf5_snr;
struct detector *det;
- IndexingMethod *indm;
- IndexingPrivate **ipriv;
+ IndexingPrivate *ipriv;
int peaks; /* Peak detection method */
float tols[4];
struct beam_params *beam;
char *hdf5_peak_path;
+ int half_pixel_shift;
float pk_inn;
float pk_mid;
float pk_out;
float ir_inn;
float ir_mid;
float ir_out;
- struct copy_hdf5_field *copyme;
+ int min_res;
+ int max_res;
+ int max_n_peaks;
+ int min_pix_count;
+ int max_pix_count;
+ int local_bg_radius;
+ int min_peaks;
+ struct imagefile_field_list *copyme;
int integrate_saturated;
int use_saturated;
int no_revalidate;
int stream_peaks;
int stream_refls;
+ int stream_nonhits;
IntegrationMethod int_meth;
IntDiag int_diag;
signed int int_diag_h;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 20:13:32 +0000