diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/ambigator.c | 33 | ||||
-rw-r--r-- | src/cell_explorer.c | 5 | ||||
-rw-r--r-- | src/compare_hkl.c | 119 | ||||
-rw-r--r-- | src/diffraction.c | 51 | ||||
-rw-r--r-- | src/dw-hdfsee.c | 903 | ||||
-rw-r--r-- | src/dw-hdfsee.h | 13 | ||||
-rw-r--r-- | src/geoptimiser.c | 1707 | ||||
-rw-r--r-- | src/hdfsee-render.c | 38 | ||||
-rw-r--r-- | src/hdfsee.c | 13 | ||||
-rw-r--r-- | src/hrs-scaling.c | 559 | ||||
-rw-r--r-- | src/im-sandbox.c | 102 | ||||
-rw-r--r-- | src/indexamajig.c | 63 | ||||
-rw-r--r-- | src/list_events.c | 195 | ||||
-rw-r--r-- | src/merge.c | 248 | ||||
-rw-r--r-- | src/merge.h (renamed from src/hrs-scaling.h) | 21 | ||||
-rw-r--r-- | src/partial_sim.c | 32 | ||||
-rw-r--r-- | src/partialator.c | 424 | ||||
-rw-r--r-- | src/pattern_sim.c | 5 | ||||
-rw-r--r-- | src/post-refinement.c | 579 | ||||
-rw-r--r-- | src/post-refinement.h | 32 | ||||
-rw-r--r-- | src/predict-refine.c | 757 | ||||
-rw-r--r-- | src/predict-refine.h | 45 | ||||
-rw-r--r-- | src/process_hkl.c | 11 | ||||
-rw-r--r-- | src/process_image.c | 191 | ||||
-rw-r--r-- | src/process_image.h | 12 | ||||
-rw-r--r-- | src/rejection.c | 191 | ||||
-rw-r--r-- | src/rejection.h | 43 | ||||
-rw-r--r-- | src/scaling-report.c | 898 | ||||
-rw-r--r-- | src/scaling-report.h | 82 |
29 files changed, 4230 insertions, 3142 deletions
diff --git a/src/ambigator.c b/src/ambigator.c index 39609e5f..2268b400 100644 --- a/src/ambigator.c +++ b/src/ambigator.c @@ -3,13 +3,13 @@ * * Resolve indexing ambiguities * - * Copyright © 2014 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2014 Wolfgang Brehm * * Authors: - * 2014 Thomas White <taw@physics.org> - * 2014 Wolfgang Brehm <wolfgang.brehm@gmail.com> + * 2014-2015 Thomas White <taw@physics.org> + * 2014 Wolfgang Brehm <wolfgang.brehm@gmail.com> * * This file is part of CrystFEL. * @@ -1049,18 +1049,6 @@ int main(int argc, char *argv[]) } - if ( argc != (optind+1) ) { - ERROR("Please provide exactly one stream filename.\n"); - return 1; - } - - infile = argv[optind++]; - st = open_stream_for_read(infile); - if ( st == NULL ) { - ERROR("Failed to open input stream '%s'\n", infile); - return 1; - } - if ( s_sym_str == NULL ) { ERROR("You must specify the input symmetry (with -y)\n"); return 1; @@ -1111,6 +1099,19 @@ int main(int argc, char *argv[]) } } + if ( argc != (optind+1) ) { + ERROR("Please provide exactly one stream filename.\n"); + return 1; + } + + infile = argv[optind++]; + st = open_stream_for_read(infile); + if ( st == NULL ) { + ERROR("Failed to open input stream '%s'\n", infile); + return 1; + } + + crystals = NULL; n_crystals = 0; max_crystals = 0; diff --git a/src/cell_explorer.c b/src/cell_explorer.c index 93f4a26f..923462bd 100644 --- a/src/cell_explorer.c +++ b/src/cell_explorer.c @@ -418,7 +418,7 @@ static gboolean draw_sig(GtkWidget *da, GdkEventExpose *event, HistoBox *b) cairo_text_extents_t ext; char label[256]; - cairo_select_font_face(cr, "Serif", CAIRO_FONT_SLANT_ITALIC, + cairo_select_font_face(cr, "Liberation Serif", CAIRO_FONT_SLANT_ITALIC, CAIRO_FONT_WEIGHT_BOLD); cairo_set_font_size(cr, height/10.0); @@ -695,7 +695,8 @@ static void scan_minmax(CellWindow *w) fprintf(stderr, "Too many indexing methods\n"); } else { IndexingMethod m = w->indms[i]; - w->unique_indms[w->n_unique_indms++] = m; + w->unique_indms[w->n_unique_indms] = m; + w->active_indms[w->n_unique_indms++] = 1; } } diff --git a/src/compare_hkl.c b/src/compare_hkl.c index eed92af0..abb84e43 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -41,6 +41,7 @@ #include <assert.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_statistics.h> +#include <gsl/gsl_fit.h> #include "version.h" #include "utils.h" @@ -661,6 +662,110 @@ static int get_bin(struct shells *s, Reflection *refl, UnitCell *cell) } +static int wilson_scale(RefList *list1, RefList *list2, UnitCell *cell) +{ + Reflection *refl1; + Reflection *refl2; + RefListIterator *iter; + int max_n = 256; + int n = 0; + double *x; + double *y; + int r; + double G, B; + double c0, c1, cov00, cov01, cov11, chisq; + + x = malloc(max_n*sizeof(double)); + y = malloc(max_n*sizeof(double)); + if ( (x==NULL) || (y==NULL) ) { + ERROR("Failed to allocate memory for scaling.\n"); + return 1; + } + + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) + { + signed int h, k, l; + double Ih1, Ih2; + double res; + + get_indices(refl1, &h, &k, &l); + res = resolution(cell, h, k, l); + + refl2 = find_refl(list2, h, k, l); + assert(refl2 != NULL); + + Ih1 = get_intensity(refl1); + Ih2 = get_intensity(refl2); + + if ( (Ih1 <= 0.0) || (Ih2 <= 0.0) ) continue; + if ( isnan(Ih1) || isinf(Ih1) ) continue; + if ( isnan(Ih2) || isinf(Ih2) ) continue; + + if ( n == max_n ) { + max_n *= 2; + x = realloc(x, max_n*sizeof(double)); + y = realloc(y, max_n*sizeof(double)); + if ( (x==NULL) || (y==NULL) ) { + ERROR("Failed to allocate memory for scaling.\n"); + return 1; + } + } + + x[n] = res*res; + y[n] = log(Ih1/Ih2); + n++; + + } + + if ( n < 2 ) { + ERROR("Not enough reflections for scaling\n"); + return 1; + } + + r = gsl_fit_linear(x, 1, y, 1, n, &c0, &c1, + &cov00, &cov01, &cov11, &chisq); + + if ( r ) { + ERROR("Scaling failed.\n"); + return 1; + } + + G = exp(c0); + B = c1/2.0; + + STATUS("Relative scale factor = %f, relative B factor = %f A^2\n", + G, B*1e20); + STATUS("A scale factor greater than 1 means that the second reflection " + "list is weaker than the first.\n"); + STATUS("A positive relative B factor means that the second reflection " + "list falls off with resolution more quickly than the first.\n"); + + free(x); + free(y); + + /* Apply the scaling factor */ + for ( refl2 = first_refl(list2, &iter); + refl2 != NULL; + refl2 = next_refl(refl2, iter) ) + { + signed int h, k, l; + double res; + double corr; + + get_indices(refl2, &h, &k, &l); + res = resolution(cell, h, k, l); + + corr = G * exp(2.0*B*res*res); + set_intensity(refl2, get_intensity(refl2)*corr); + set_esd_intensity(refl2, get_esd_intensity(refl2)*corr); + + } + return 0; +} + + static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, double rmin, double rmax, enum fom fom, int config_unity, int nshells, const char *filename, @@ -671,7 +776,6 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, Reflection *refl1; RefListIterator *iter; FILE *fh; - double scale; struct fom_context *fctx; struct shells *shells; const char *t1, *t2; @@ -684,10 +788,9 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, return; } - if ( config_unity ) { - scale = 1.0; - } else { - scale = stat_scale_intensity(list1, list2); + if ( !config_unity && wilson_scale(list1, list2, cell) ) { + ERROR("Error with scaling.\n"); + return; } /* Calculate the bins */ @@ -726,9 +829,9 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, } i1 = get_intensity(refl1); - i2 = scale * get_intensity(refl2); + i2 = get_intensity(refl2); sig1 = get_esd_intensity(refl1); - sig2 = scale * get_esd_intensity(refl2); + sig2 = get_esd_intensity(refl2); if ( (fom == FOM_CCANO) || (fom == FOM_CRDANO) || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) ) @@ -754,7 +857,7 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, assert(refl2_bij != NULL); i1bij = get_intensity(refl1_bij); - i2bij = scale * get_intensity(refl2_bij); + i2bij = get_intensity(refl2_bij); } else { diff --git a/src/diffraction.c b/src/diffraction.c index 1fb63ea3..5a936809 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -423,6 +423,40 @@ static void diffraction_at_k(struct image *image, const double *intensities, } +static void normalise_sampled_spectrum(struct sample *spec, int n, int nsamp) +{ + int i; + double total_w = 0.0; + double samp_w = 0.0; + + for ( i=0; i<n; i++ ) { + total_w += spec[i].weight; + } + STATUS("%i energies in spectrum, %i samples requested.\n", n, nsamp); + + for ( i=0; i<nsamp; i++ ) { + samp_w += spec[i].weight; + } + + if ( samp_w < 0.8 * total_w ) { + ERROR("WARNING: Only %.1f %% of the calculated spectrum " + "intensity was sampled.\n", 100.0 * samp_w / total_w); + ERROR("I've increased the weighting of the sampled points to " + "keep the final intensity constant, but the spectrum " + "might be very far from accurate. Consider increasing " + "the number of spectrum samples.\n"); + } else { + STATUS("%.1f %% of calculated spectrum intensity sampled.\n", + 100.0 * samp_w / total_w); + STATUS("Normalised to keep total intensity constant.\n"); + } + + for ( i=0; i<n; i++ ) { + spec[i].weight /= samp_w; + } +} + + static int compare_samples(const void *a, const void *b) { struct sample *sample1 = (struct sample *)a; @@ -548,6 +582,7 @@ struct sample *generate_tophat(struct image *image) } image->spectrum_size = image->nsamples; + /* No need to call normalise_sampled_spectrum() in this case */ return spectrum; } @@ -556,7 +591,6 @@ struct sample *generate_tophat(struct image *image) struct sample *generate_SASE(struct image *image, gsl_rng *rng) { struct sample *spectrum; - int i; const int spec_size = 1024; double eV_cen; /* Central photon energy for this spectrum */ const double jitter_sigma_eV = 8.0; @@ -581,21 +615,13 @@ struct sample *generate_SASE(struct image *image, gsl_rng *rng) /* Add SASE-type noise to Gaussian spectrum */ add_sase_noise(spectrum, spec_size, rng); - /* Normalise intensity (before taking restricted number of samples) */ - double total_weight = 0.0; - for ( i=0; i<spec_size; i++ ) { - total_weight += spectrum[i].weight; - } - for ( i=0; i<spec_size; i++ ) { - spectrum[i].weight /= total_weight; - } - /* Sort samples in spectrum by weight. Diffraction calculation will * take the requested number, starting from the brightest */ qsort(spectrum, spec_size, sizeof(struct sample), compare_samples); - image->spectrum_size = spec_size; + normalise_sampled_spectrum(spectrum, spec_size, image->nsamples); + image->spectrum_size = spec_size; return spectrum; } @@ -662,8 +688,9 @@ struct sample *generate_twocolour(struct image *image) * take the requested number, starting from the brightest */ qsort(spectrum, spec_size, sizeof(struct sample), compare_samples); - image->spectrum_size = spec_size; + normalise_sampled_spectrum(spectrum, spec_size, image->nsamples); + image->spectrum_size = spec_size; return spectrum; } diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c index e92dcd1e..a899acc5 100644 --- a/src/dw-hdfsee.c +++ b/src/dw-hdfsee.c @@ -144,7 +144,7 @@ static void draw_panel_rectangle(cairo_t *cr, cairo_matrix_t *basic_m, } -int render_adsc_uint16(DisplayWindow *dw, const char *filename) +static int render_adsc_uint16(DisplayWindow *dw, const char *filename) { int x, y, fs, ss; double dfs, dss; @@ -853,9 +853,9 @@ static gint displaywindow_set_boostint(GtkWidget *widget, DisplayWindow *dw) return 0; } - if ( dw->hdfile == NULL ) { - return 0; - } + if ( dw->hdfile == NULL ) { + return 0; + } bd = malloc(sizeof(BoostIntDialog)); if ( bd == NULL ) return 0; @@ -904,6 +904,207 @@ static gint displaywindow_set_boostint(GtkWidget *widget, DisplayWindow *dw) } +static void do_filters(DisplayWindow *dw) +{ + if ( dw->median_filter > 0 ) { + filter_median(dw->image, dw->median_filter); + } + + if ( dw->noisefilter ) { + filter_noise(dw->image); + } +} + + +static gint displaywindow_newevent(DisplayWindow *dw, int new_event) +{ + int fail; + int i; + char title[1024]; + char *bn; + + if ( dw->not_ready_yet ) return 0; + + float *old_data = dw->image->data; + uint16_t *old_flags = dw->image->flags; + 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); + if ( fail ) { + ERROR("Couldn't load image"); + dw->image->data = old_data; + dw->image->flags = old_flags; + dw->image->dp = old_dp; + dw->image->bad = old_bad; + return 1; + } + + dw->curr_event = new_event; + + do_filters(dw); + displaywindow_update(dw); + + bn = safe_basename(dw->image->filename); + sprintf(title, "%s - event: %s - hdfsee", bn, + get_event_string(dw->ev_list->events[dw->curr_event])); + gtk_window_set_title(GTK_WINDOW(dw->window), title); + free(bn); + + for (i = 0; i < dw->image->det->n_panels; i++) { + free(old_dp[i]); + free(old_bad[i]); + } + free(old_data); + free(old_flags); + free(old_dp); + free(old_bad); + return 0; +} + + +static gint displaywindow_randomevent(GtkWidget *widget, DisplayWindow *dw) +{ + int rand_event; + + if ( dw->not_ready_yet ) return 0; + + rand_event = gsl_rng_uniform_int(&dw->rng, dw->ev_list->num_events); + + return displaywindow_newevent(dw, rand_event); +} + + +static gint displaywindow_set_newevent_response(GtkWidget *widget, + gint response, + DisplayWindow *dw) +{ + int ei; + int matched_event; + int done = 1; + + if ( response == GTK_RESPONSE_OK ) { + + const char *sevent; + + + + sevent = gtk_entry_get_text( + GTK_ENTRY(dw->event_dialog->entry)); + + matched_event = -1; + + for ( ei=0; ei<dw->ev_list->num_events; ei++ ) { + + char *ei_ev_string; + + ei_ev_string = get_event_string(dw->ev_list->events[ei]); + if ( strcmp(ei_ev_string, sevent) == 0 ) { + matched_event = ei; + break; + } + } + + if ( matched_event == -1 ) { + displaywindow_error(dw, "Cannot find event.\n"); + done = 0; + } else { + displaywindow_newevent(dw, matched_event); + displaywindow_update(dw); + } + } + + if ( done ) { + gtk_widget_destroy(dw->event_dialog->window); + } + + return 0; +} + + +static gint displaywindow_set_newevent_destroy(GtkWidget *widget, + DisplayWindow *dw) +{ + free(dw->event_dialog); + dw->event_dialog = NULL; + return 0; +} + + +static gint displaywindow_set_newevent_response_ac(GtkWidget *widget, + DisplayWindow *dw) +{ + return displaywindow_set_newevent_response(widget, GTK_RESPONSE_OK, dw); +} + + +/* Create a window to ask the user for a new intensity boost factor */ +static gint displaywindow_set_newevent(GtkWidget *widget, DisplayWindow *dw) +{ + EventDialog *ed; + GtkWidget *vbox; + GtkWidget *hbox; + GtkWidget *table; + GtkWidget *label; + char tmp[1024]; + + if ( dw->event_dialog != NULL ) { + return 0; + } + + if ( dw->hdfile == NULL ) { + return 0; + } + + ed = malloc(sizeof(EventDialog)); + if ( ed == NULL ) return 0; + dw->event_dialog = ed; + + ed->window = gtk_dialog_new_with_buttons("Go to event", + GTK_WINDOW(dw->window), + GTK_DIALOG_DESTROY_WITH_PARENT, + GTK_STOCK_CANCEL, GTK_RESPONSE_CLOSE, + GTK_STOCK_OK, GTK_RESPONSE_OK, NULL); + + vbox = gtk_vbox_new(FALSE, 0); + hbox = gtk_hbox_new(TRUE, 0); + gtk_box_pack_start(GTK_BOX(GTK_DIALOG(ed->window)->vbox), + GTK_WIDGET(hbox), FALSE, FALSE, 7); + gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5); + + table = gtk_table_new(3, 2, FALSE); + gtk_table_set_row_spacings(GTK_TABLE(table), 5); + gtk_table_set_col_spacings(GTK_TABLE(table), 5); + gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(table), FALSE, FALSE, 0); + + label = gtk_label_new("Event ID:"); + gtk_misc_set_alignment(GTK_MISC(label), 1, 0.5); + gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(label), + 1, 2, 3, 4); + + ed->entry = gtk_entry_new(); + snprintf(tmp, 1023, "%s", + get_event_string(dw->ev_list->events[dw->curr_event])); + gtk_entry_set_text(GTK_ENTRY(ed->entry), tmp); + gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(ed->entry), + 2, 3, 3, 4); + + g_signal_connect(G_OBJECT(ed->entry), "activate", + G_CALLBACK(displaywindow_set_newevent_response_ac), + dw); + g_signal_connect(G_OBJECT(ed->window), "response", + G_CALLBACK(displaywindow_set_newevent_response), dw); + g_signal_connect(G_OBJECT(ed->window), "destroy", + G_CALLBACK(displaywindow_set_newevent_destroy), dw); + gtk_window_set_resizable(GTK_WINDOW(ed->window), FALSE); + gtk_widget_show_all(ed->window); + gtk_widget_grab_focus(GTK_WIDGET(ed->entry)); + + return 0; +} + + static gint displaywindow_set_ringradius_response(GtkWidget *widget, gint response, DisplayWindow *dw) @@ -1140,6 +1341,23 @@ static void load_features_from_file(struct image *image, const char *filename) image_add_feature(image->features, fs, ss, image, 1.0, "peak"); + } else if ( r == 2 ) { + + p = find_orig_panel(image->det, fs, ss); + + if ( p == NULL ) { + ERROR("Unable to find panel " + "(no geometry file given?)\n"); + } else { + + /* Convert coordinates to match rearranged + * panels in memory */ + fs = fs - p->orig_min_fs + p->min_fs; + ss = ss - p->orig_min_ss + p->min_ss; + } + image_add_feature(image->features, fs, ss, image, 1.0, + "peak"); + } } while ( rval != NULL ); @@ -1191,13 +1409,15 @@ static gint displaywindow_about(GtkWidget *widget, DisplayWindow *dw) gtk_about_dialog_set_name(GTK_ABOUT_DIALOG(window), "hdfsee"); gtk_about_dialog_set_version(GTK_ABOUT_DIALOG(window), PACKAGE_VERSION); gtk_about_dialog_set_copyright(GTK_ABOUT_DIALOG(window), - "(c) 2012-2013 Thomas White <taw@physics.org> and others"); + "© 2012-2015 Deutsches Elektronen-Synchrotron DESY," + " a research centre of the Helmholtz Association."); gtk_about_dialog_set_comments(GTK_ABOUT_DIALOG(window), - "Quick viewer for HDF files"); + "Quick viewer for HDF files"); gtk_about_dialog_set_license(GTK_ABOUT_DIALOG(window), - "(c) 2012-2013 Thomas White <taw@physics.org>\n"); + "© 2012-2015 Deutsches Elektronen-Synchrotron DESY," + " a research centre of the Helmholtz Association."); gtk_about_dialog_set_website(GTK_ABOUT_DIALOG(window), - "http://www.desy.de/~twhite/crystfel"); + "http://www.desy.de/~twhite/crystfel"); gtk_about_dialog_set_authors(GTK_ABOUT_DIALOG(window), authors); g_signal_connect(window, "response", G_CALLBACK(gtk_widget_destroy), @@ -1209,7 +1429,7 @@ static gint displaywindow_about(GtkWidget *widget, DisplayWindow *dw) } -static int save_geometry_file(DisplayWindow *dw) +static int save_geometry_file(GtkWidget *widget, DisplayWindow *dw) { GtkWidget *d; gchar *output_filename; @@ -1227,9 +1447,11 @@ static int save_geometry_file(DisplayWindow *dw) w = write_detector_geometry_2(dw->geom_filename, output_filename, dw->image->det, "Manually optimized with " "hdfsee", 0); - if ( w != 0 ) { - ERROR("Error saving geometry!\n"); + if ( w != 0 && w!=2 ) { + displaywindow_error(dw, + "Unable to save the detector geometry."); } + gtk_widget_destroy(d); g_free(output_filename); return w; @@ -1256,6 +1478,54 @@ static gint displaywindow_peak_overlay(GtkWidget *widget, DisplayWindow *dw) } +static void set_calibration_menu_sensitivity(DisplayWindow *dw, int val) { + + GtkAction * a; + + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/calibration"); + gtk_action_set_sensitive(GTK_ACTION(a), val); + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/calibration/calibrationprevious"); + gtk_action_set_sensitive(GTK_ACTION(a), val); + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/calibration/calibrationnext"); + gtk_action_set_sensitive(GTK_ACTION(a), val); + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/calibration/switchcalibmode"); + gtk_action_set_sensitive(GTK_ACTION(a), val); + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/calibration/focus"); + gtk_action_set_sensitive(GTK_ACTION(a), val); + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/calibration/savegeometry"); + gtk_action_set_sensitive(GTK_ACTION(a), val); +} + + +static void set_events_menu_sensitivity(DisplayWindow *dw, int val) { + + GtkAction * a; + + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/events"); + gtk_action_set_sensitive(GTK_ACTION(a), val); + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/events/eventprevious"); + gtk_action_set_sensitive(GTK_ACTION(a), val); + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/events/eventnext"); + gtk_action_set_sensitive(GTK_ACTION(a), val); + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/events/gotoevent"); + gtk_action_set_sensitive(GTK_ACTION(a), val); + a = gtk_ui_manager_get_action(dw->ui, + "/ui/displaywindow/events/randomevent"); + gtk_action_set_sensitive(GTK_ACTION(a), val); + +} + + static gint displaywindow_set_calibmode(GtkWidget *d, DisplayWindow *dw) { GtkWidget *w, *vbox; @@ -1263,7 +1533,7 @@ static gint displaywindow_set_calibmode(GtkWidget *d, DisplayWindow *dw) w = gtk_ui_manager_get_widget(dw->ui, "/ui/displaywindow/tools/calibmode"); - if ( dw->image->det == dw->simple_geom ) { + if ( dw->simple ) { gtk_check_menu_item_set_state(GTK_CHECK_MENU_ITEM(w), 0); return 0; } @@ -1303,6 +1573,8 @@ static gint displaywindow_set_calibmode(GtkWidget *d, DisplayWindow *dw) dw->calib_mode_curr_p = dw->calib_mode_curr_rg->panels[0]; } + set_calibration_menu_sensitivity(dw, TRUE); + dw->calib_mode = CALIBMODE_PANELS; dw->statusbar = gtk_statusbar_new(); @@ -1318,6 +1590,8 @@ static gint displaywindow_set_calibmode(GtkWidget *d, DisplayWindow *dw) } else { + set_calibration_menu_sensitivity(dw, FALSE); + dw->calib_mode = CALIBMODE_NONE; gtk_widget_destroy(dw->statusbar); dw->statusbar = NULL; @@ -1694,6 +1968,207 @@ static gint displaywindow_setscale(GtkWidget *widget, GtkRadioAction *action, return 0; } +static int curr_rg_pointer_index(DisplayWindow *dw) +{ + int r; + + for ( r=0; r<dw->rg_coll->n_rigid_groups; ++r) { + if ( dw->rg_coll->rigid_groups[r] == dw->calib_mode_curr_rg ) { + return r; + } + } + + /* Never reached (we hope) */ + return 999; +} + +static int curr_p_pointer_index(DisplayWindow *dw) +{ + int p; + + for ( p=0; p<dw->image->det->n_panels; ++p) { + if ( &dw->image->det->panels[p] == dw->calib_mode_curr_p ) { + return p; + } + } + + /* Never reached (we hope) */ + return 999; +} + + +static void select_next_group(DisplayWindow *dw, int num_rg) +{ + if ( dw->calib_mode_curr_rg == dw->rg_coll->rigid_groups[num_rg-1] ) { + dw->calib_mode_curr_rg = dw->rg_coll->rigid_groups[0]; + } else { + dw->calib_mode_curr_rg = + dw->rg_coll->rigid_groups[curr_rg_pointer_index(dw)+1]; + } +} + + +static void select_prev_group(DisplayWindow *dw, int num_rg) +{ + if ( dw->calib_mode_curr_rg == dw->rg_coll->rigid_groups[0] ) { + dw->calib_mode_curr_rg = dw->rg_coll->rigid_groups[num_rg-1]; + } else { + dw->calib_mode_curr_rg = + dw->rg_coll->rigid_groups[curr_rg_pointer_index(dw)-1]; + } +} + + +static void select_next_panel(DisplayWindow *dw, int num_p) +{ + if ( dw->calib_mode_curr_p == &dw->image->det->panels[num_p-1] ) { + dw->calib_mode_curr_p = &dw->image->det->panels[0]; + } else { + dw->calib_mode_curr_p = + &dw->image->det->panels[curr_p_pointer_index(dw)+1]; + } +} + + +static void select_prev_panel(DisplayWindow *dw, int num_p) +{ + if ( dw->calib_mode_curr_p == &dw->image->det->panels[0] ) { + dw->calib_mode_curr_p = &dw->image->det->panels[num_p-1]; + } else { + dw->calib_mode_curr_p = + &dw->image->det->panels[curr_p_pointer_index(dw)-1]; + } +} + + +static void toggle_calibmode_groupmode(GtkWidget *widget, DisplayWindow *dw) +{ + struct rigid_group *rg; + struct detector *det = dw->image->det; + + switch ( dw->calib_mode ) { + + case CALIBMODE_NONE: + break; + + case CALIBMODE_PANELS: + if ( det->n_rigid_groups != det->n_panels ) { + /* Only change if there are any rigid groups defined */ + dw->calib_mode = CALIBMODE_GROUPS; + rg = find_corresponding_rigid_group(dw, + dw->calib_mode_curr_p); + if ( rg == NULL) { + dw->calib_mode = CALIBMODE_ALL; + } else { + dw->calib_mode_curr_rg = rg; + } + + } else { + /* ...otherwise skip to ALL mode */ + dw->calib_mode = CALIBMODE_ALL; + } + break; + + case CALIBMODE_GROUPS: + dw->calib_mode = CALIBMODE_ALL; + break; + + case CALIBMODE_ALL: + dw->calib_mode = CALIBMODE_PANELS; + dw->calib_mode_curr_p = dw->calib_mode_curr_rg->panels[0]; + break; + + } + redraw_window(dw); +} + + +static void toggle_calibmode_focus(GtkWidget *widget, DisplayWindow *dw) +{ + dw->calib_mode_show_focus = 1 - dw->calib_mode_show_focus; + redraw_window(dw); +} + + +static void calibmode_next(GtkWidget *widget, DisplayWindow *dw) +{ + int n; + + switch ( dw->calib_mode ) { + + case CALIBMODE_NONE: + break; + + case CALIBMODE_PANELS: + n = dw->image->det->n_panels; + select_next_panel(dw, n); + break; + + case CALIBMODE_GROUPS: + n = dw->image->det->n_rigid_groups; + select_next_group(dw, n); + break; + + case CALIBMODE_ALL: + break; + + } + redraw_window(dw); +} + + +static void calibmode_prev(GtkWidget *widget, DisplayWindow *dw) +{ + int n; + + switch ( dw->calib_mode ) { + + case CALIBMODE_NONE: + break; + + case CALIBMODE_PANELS: + n = dw->image->det->n_panels; + select_prev_panel(dw, n); + break; + + case CALIBMODE_GROUPS: + n = dw->image->det->n_rigid_groups; + select_prev_group(dw, n); + break; + + case CALIBMODE_ALL: + break; + + } + redraw_window(dw); +} + + +static void event_next(GtkWidget *widget, DisplayWindow *dw) +{ + int new_event; + + if ( dw->curr_event == dw->ev_list->num_events-1 ) { + new_event = 0; + } else { + new_event = dw->curr_event+1; + } + displaywindow_newevent(dw, new_event); +} + + +static void event_prev(GtkWidget *widget, DisplayWindow *dw) +{ + int new_event; + + if ( dw->curr_event == 0 ) { + new_event = dw->ev_list->num_events-1; + } else { + new_event = dw->curr_event-1; + } + displaywindow_newevent(dw, new_event); +} + static void displaywindow_addmenubar(DisplayWindow *dw, GtkWidget *vbox, int colscale) @@ -1722,7 +2197,27 @@ static void displaywindow_addmenubar(DisplayWindow *dw, GtkWidget *vbox, { "PeaksAction", NULL, "Load Feature List...", NULL, NULL, G_CALLBACK(displaywindow_peak_overlay) }, - { "EventsAction", NULL, "_Events", NULL, NULL, NULL }, + { "CalibrationAction", NULL, "_Calibration", NULL, NULL, NULL }, + { "CalibPreviousAction", NULL, "Previous Item", "minus", NULL, + G_CALLBACK(calibmode_prev) }, + { "CalibNextAction", NULL, "Next Item", "plus", NULL, + G_CALLBACK(calibmode_next) }, + { "SwitchCalibModeAction", NULL, "Toggle Panel/Group/All", "g", + NULL, G_CALLBACK(toggle_calibmode_groupmode) }, + { "ToggleFocusAction", NULL, "Toggle Focus Rectangle", "i", + NULL, G_CALLBACK(toggle_calibmode_focus) }, + { "SaveGeometryAction", NULL, "Save Geometry", "s", NULL, + G_CALLBACK(save_geometry_file) }, + + { "EventsAction", NULL, "_Event", NULL, NULL, NULL }, + { "EventPreviousAction", NULL, "Previous", "p", NULL, + G_CALLBACK(event_prev) }, + { "EventNextAction", NULL, "Next", "n", NULL, + G_CALLBACK(event_next) }, + { "GotoEventAction", NULL, "Go To Event", "e", NULL, + G_CALLBACK(displaywindow_set_newevent) }, + { "RandomEventAction", NULL, "Go To Random Event", "r", NULL, + G_CALLBACK(displaywindow_randomevent) }, { "HelpAction", NULL, "_Help", NULL, NULL, NULL }, { "AboutAction", GTK_STOCK_ABOUT, "_About hdfsee...", @@ -1779,17 +2274,6 @@ static void displaywindow_addmenubar(DisplayWindow *dw, GtkWidget *vbox, } -static void do_filters(DisplayWindow *dw) -{ - if ( dw->median_filter > 0 ) { - filter_median(dw->image, dw->median_filter); - } - - if ( dw->noisefilter ) { - filter_noise(dw->image); - } -} - struct newhdf { DisplayWindow *dw; GtkWidget *widget; @@ -1806,15 +2290,16 @@ static gint displaywindow_newhdf(GtkMenuItem *item, struct newhdf *nh) a = gtk_check_menu_item_get_active(GTK_CHECK_MENU_ITEM(nh->widget)); if ( !a ) return 0; + /* hdf5_read() will create a new simple geom, so get rid of the old + * 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); if ( fail ) { ERROR("Couldn't load image"); return 1; } - nh->dw->simple_geom = simple_geometry(nh->dw->image); - nh->dw->image->det = nh->dw->simple_geom; - do_filters(nh->dw); displaywindow_update(nh->dw); return 0; @@ -1958,102 +2443,6 @@ static int displaywindow_update_menus(DisplayWindow *dw, const char *selectme) } -static gint displaywindow_newevent(GtkMenuItem *item, struct newev *ne) -{ - gboolean a; - int fail; - int i; - - if ( ne->dw->not_ready_yet ) return 0; - - a = gtk_check_menu_item_get_active(GTK_CHECK_MENU_ITEM(ne->widget)); - if ( !a ) return 0; - - float *old_data = ne->dw->image->data; - uint16_t *old_flags = ne->dw->image->flags; - float **old_dp = ne->dw->image->dp; - int **old_bad = ne->dw->image->bad; - - fail = hdf5_read2(ne->dw->hdfile, ne->dw->image, - ne->dw->ev_list->events[ne->new_ev], 0); - if ( fail ) { - ERROR("Couldn't load image"); - return 1; - } - - ne->dw->curr_event = ne->new_ev; - - do_filters(ne->dw); - displaywindow_update(ne->dw); - for (i = 0; i < ne->dw->image->det->n_panels; i++) { - free(old_dp[i]); - free(old_bad[i]); - } - free(old_data); - free(old_flags); - free(old_dp); - free(old_bad); - return 0; -} - - -static int displaywindow_update_event_menu(DisplayWindow *dw, - struct event_list *ev_list, - int curr_event) -{ - - int ei; - GtkWidget *w; - GtkWidget *ww; - GSList *grp = NULL; - - w = gtk_ui_manager_get_widget(dw->ui, "/ui/displaywindow/events"); - ww = gtk_menu_new(); - - for ( ei=0; ei<ev_list->num_events; ei++ ) { - - GtkWidget *www; - struct newev *ne; - char *ev_string; - - ev_string = get_event_string(ev_list->events[ei]); - www = gtk_radio_menu_item_new_with_label(grp, ev_string); - free(ev_string); - - ne = malloc(sizeof(struct newev)); - if ( ne != NULL ) { - - ne->widget = www; - ne->dw = dw; - ne->new_ev = ei; - - g_signal_connect(G_OBJECT(www), "toggled", - G_CALLBACK(displaywindow_newevent), ne); - - } - - if ( ei == dw->curr_event ) { - gtk_check_menu_item_set_active(GTK_CHECK_MENU_ITEM(www), - TRUE); - } else { - gtk_check_menu_item_set_active(GTK_CHECK_MENU_ITEM(www), - FALSE); - } - - gtk_menu_shell_append(GTK_MENU_SHELL(ww), www); - - grp = gtk_radio_menu_item_get_group(GTK_RADIO_MENU_ITEM(www)); - gtk_widget_show_all(www); - - } - - gtk_menu_item_set_submenu(GTK_MENU_ITEM(w), ww); - - return 0; -} - - - static gint displaywindow_release(GtkWidget *widget, GdkEventButton *event, DisplayWindow *dw) { @@ -2151,174 +2540,6 @@ static gint displaywindow_press(GtkWidget *widget, GdkEventButton *event, return 0; } - -static int curr_rg_pointer_index(DisplayWindow *dw) -{ - int r; - - for ( r=0; r<dw->rg_coll->n_rigid_groups; ++r) { - if ( dw->rg_coll->rigid_groups[r] == dw->calib_mode_curr_rg ) { - return r; - } - } - - /* Never reached (we hope) */ - return 999; -} - - -static int curr_p_pointer_index(DisplayWindow *dw) -{ - int p; - - for ( p=0; p<dw->image->det->n_panels; ++p) { - if ( &dw->image->det->panels[p] == dw->calib_mode_curr_p ) { - return p; - } - } - - /* Never reached (we hope) */ - return 999; -} - - -static void select_next_group(DisplayWindow *dw, int num_rg) -{ - if ( dw->calib_mode_curr_rg == dw->rg_coll->rigid_groups[num_rg-1] ) { - dw->calib_mode_curr_rg = dw->rg_coll->rigid_groups[0]; - } else { - dw->calib_mode_curr_rg = - dw->rg_coll->rigid_groups[curr_rg_pointer_index(dw)+1]; - } -} - - -static void select_prev_group(DisplayWindow *dw, int num_rg) -{ - if ( dw->calib_mode_curr_rg == dw->rg_coll->rigid_groups[0] ) { - dw->calib_mode_curr_rg = dw->rg_coll->rigid_groups[num_rg-1]; - } else { - dw->calib_mode_curr_rg = - dw->rg_coll->rigid_groups[curr_rg_pointer_index(dw)-1]; - } -} - - -static void select_next_panel(DisplayWindow *dw, int num_p) -{ - if ( dw->calib_mode_curr_p == &dw->image->det->panels[num_p-1] ) { - dw->calib_mode_curr_p = &dw->image->det->panels[0]; - } else { - dw->calib_mode_curr_p = - &dw->image->det->panels[curr_p_pointer_index(dw)+1]; - } -} - - -static void select_prev_panel(DisplayWindow *dw, int num_p) -{ - if ( dw->calib_mode_curr_p == &dw->image->det->panels[0] ) { - dw->calib_mode_curr_p = &dw->image->det->panels[num_p-1]; - } else { - dw->calib_mode_curr_p = - &dw->image->det->panels[curr_p_pointer_index(dw)-1]; - } -} - - -static void toggle_calibmode_groupmode(DisplayWindow *dw) -{ - struct rigid_group *rg; - struct detector *det = dw->image->det; - - switch ( dw->calib_mode ) { - - case CALIBMODE_NONE: - break; - - case CALIBMODE_PANELS: - if ( det->n_rigid_groups != det->n_panels ) { - /* Only change if there are any rigid groups defined */ - dw->calib_mode = CALIBMODE_GROUPS; - rg = find_corresponding_rigid_group(dw, - dw->calib_mode_curr_p); - if ( rg == NULL) { - dw->calib_mode = CALIBMODE_ALL; - } else { - dw->calib_mode_curr_rg = rg; - } - - } else { - /* ...otherwise skip to ALL mode */ - dw->calib_mode = CALIBMODE_ALL; - } - break; - - case CALIBMODE_GROUPS: - dw->calib_mode = CALIBMODE_ALL; - break; - - case CALIBMODE_ALL: - dw->calib_mode = CALIBMODE_PANELS; - dw->calib_mode_curr_p = dw->calib_mode_curr_rg->panels[0]; - break; - - } -} - - -static void calibmode_next(DisplayWindow *dw) -{ - int n; - - switch ( dw->calib_mode ) { - - case CALIBMODE_NONE: - break; - - case CALIBMODE_PANELS: - n = dw->image->det->n_panels; - select_next_panel(dw, n); - break; - - case CALIBMODE_GROUPS: - n = dw->image->det->n_rigid_groups; - select_next_group(dw, n); - break; - - case CALIBMODE_ALL: - break; - - } -} - - -static void calibmode_prev(DisplayWindow *dw) -{ - int n; - - switch ( dw->calib_mode ) { - - case CALIBMODE_NONE: - break; - - case CALIBMODE_PANELS: - n = dw->image->det->n_panels; - select_prev_panel(dw, n); - break; - - case CALIBMODE_GROUPS: - n = dw->image->det->n_rigid_groups; - select_prev_group(dw, n); - break; - - case CALIBMODE_ALL: - break; - - } -} - - static void calibmode_up(DisplayWindow *dw) { int pi; @@ -2438,7 +2659,6 @@ static void calibmode_right(DisplayWindow *dw) static gint displaywindow_keypress(GtkWidget *widget, GdkEventKey *event, DisplayWindow *dw) { - int s; if ( !dw->calib_mode ) { return 0; @@ -2470,41 +2690,40 @@ static gint displaywindow_keypress(GtkWidget *widget, GdkEventKey *event, redraw_window(dw); break; - case GDK_plus: case GDK_KP_Add: - calibmode_next(dw); - redraw_window(dw); + calibmode_next(NULL, dw); break; - case GDK_minus: case GDK_KP_Subtract: - calibmode_prev(dw); - redraw_window(dw); + calibmode_prev(NULL, dw); break; + } - case GDK_f: - dw->calib_mode_show_focus = 1 - dw->calib_mode_show_focus; - redraw_window(dw); - break; + return 0; +} - case GDK_g: - toggle_calibmode_groupmode(dw); - redraw_window(dw); - break; - case GDK_s: - s = save_geometry_file(dw); - if ( s != 0 ) { - if ( s != 2 ) { - displaywindow_error(dw, - "Unable to save the detector geometry."); - } - } - break; +static void impose_twod_geometry(DisplayWindow *dw, const char *twod_element) +{ + + int i; + + for ( i=0; i<dw->image->det->n_panels; i++ ) { + + struct panel *p; + + p = &dw->image->det->panels[i]; + if ( p->data != NULL ) free(p->data); + p->data = strdup(twod_element); + + if ( p->dim_structure ) free_dim_structure(p->dim_structure); + p->dim_structure = default_dim_structure(); } - return 0; + dw->image->det->path_dim = 0; + dw->image->det->dim_dim = 0; + } @@ -2521,11 +2740,10 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, int median_filter) { DisplayWindow *dw; - char *title; GtkWidget *vbox; int check; GtkWidget *w; - GtkWidget *ww; + char title[1024]; dw = calloc(1, sizeof(DisplayWindow)); if ( dw == NULL ) return NULL; @@ -2537,7 +2755,7 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, dw->boostint = 1; dw->motion_callback = 0; dw->numbers_window = NULL; - dw->simple_geom = NULL; + dw->simple = 0; dw->image = NULL; dw->show_rings = show_rings; dw->show_peaks = 0; @@ -2559,6 +2777,14 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, dw->statusbar = NULL; dw->multi_event = 0; dw->ev_list = NULL; + dw->rng = *gsl_rng_alloc(gsl_rng_mt19937); + FILE *fh; + unsigned long int seed; + fh = fopen("/dev/urandom", "r"); + fread(&seed, sizeof(seed), 1, fh); + fclose(fh); + gsl_rng_set(&dw->rng, seed); + if ( geom_filename != NULL ) { dw->geom_filename = strdup(geom_filename); } else { @@ -2583,6 +2809,11 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, return NULL; } + if ( dw->image->det != NULL && element != NULL ) { + impose_twod_geometry(dw, element); + dw->multi_event = 0; + } + if ( dw->image->det != NULL && ( dw->image->det->path_dim != 0 || dw->image->det->dim_dim != 0 )) { @@ -2633,6 +2864,7 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, } else { check = hdf5_read(dw->hdfile, dw->image, element, 0); + dw->simple = 1; } if ( check ) { ERROR("Couldn't load file\n"); @@ -2643,11 +2875,6 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, dw->image->filename = strdup(filename); - if ( dw->image->det == NULL ) { - dw->simple_geom = simple_geometry(dw->image); - dw->image->det = dw->simple_geom; - } - /* Filters need geometry */ do_filters(dw); @@ -2659,11 +2886,15 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, dw->window = gtk_window_new(GTK_WINDOW_TOPLEVEL); char *bn = safe_basename(filename); - title = malloc(strlen(bn)+14); - sprintf(title, "%s - hdfsee", bn); + if ( dw->multi_event == 0 ) { + sprintf(title, "%s - hdfsee", bn); + } else { + sprintf(title, "%s - event: %s - hdfsee", bn, + get_event_string(dw->ev_list->events[dw->curr_event])); + } + free(bn); gtk_window_set_title(GTK_WINDOW(dw->window), title); - free(title); g_signal_connect(G_OBJECT(dw->window), "destroy", G_CALLBACK(displaywindow_closed), dw); @@ -2690,17 +2921,16 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, w = gtk_ui_manager_get_widget(dw->ui, "/ui/displaywindow/view/images"); - if ( dw->image->det != dw->simple_geom ) { + if ( !dw->simple ) { gtk_widget_set_sensitive(GTK_WIDGET(w), FALSE); } - ww = gtk_ui_manager_get_widget(dw->ui, - "/ui/displaywindow/events"); - - if ( dw->image->det == dw->simple_geom || dw->multi_event == 0) { - gtk_widget_set_sensitive(GTK_WIDGET(ww), FALSE); + if ( dw->simple || dw->multi_event == 0) { + set_events_menu_sensitivity(dw, FALSE); } + set_calibration_menu_sensitivity(dw, FALSE); + displaywindow_update(dw); gtk_widget_add_events(GTK_WIDGET(dw->drawingarea), @@ -2719,13 +2949,8 @@ DisplayWindow *displaywindow_open(char *filename, char *geom_filename, g_signal_connect(GTK_OBJECT(dw->drawingarea), "key-press-event", G_CALLBACK(displaywindow_keypress), dw); - if ( dw->image->det == dw->simple_geom ) { + if ( dw->simple ) { displaywindow_update_menus(dw, element); - } else { - if ( dw->multi_event ) { - displaywindow_update_event_menu(dw, dw->ev_list, - dw->curr_event); - } } dw->not_ready_yet = 0; diff --git a/src/dw-hdfsee.h b/src/dw-hdfsee.h index 81108213..85c82ac7 100644 --- a/src/dw-hdfsee.h +++ b/src/dw-hdfsee.h @@ -59,6 +59,12 @@ typedef struct { } RingRadiusDialog; +typedef struct { + GtkWidget *window; + GtkWidget *entry; +} EventDialog; + + struct numberswindow { GtkWidget *window; GtkWidget *labels[17*17]; @@ -83,6 +89,9 @@ typedef struct { GtkWidget *scrollarea; GtkUIManager *ui; GtkActionGroup *action_group; + GtkActionGroup *calibration_action_group; + GtkActionGroup *events_action_group; + int n_pixbufs; GdkPixbuf **pixbufs; gulong motion_callback; @@ -90,7 +99,7 @@ typedef struct { int not_ready_yet; - struct detector *simple_geom; + int simple; struct hdfile *hdfile; struct image *image; @@ -102,6 +111,7 @@ typedef struct { BinningDialog *binning_dialog; BoostIntDialog *boostint_dialog; RingRadiusDialog *ringradius_dialog; + EventDialog *event_dialog; struct numberswindow *numbers_window; int width; @@ -136,6 +146,7 @@ typedef struct { struct event_list *ev_list; int curr_event; + gsl_rng rng; } DisplayWindow; diff --git a/src/geoptimiser.c b/src/geoptimiser.c index 002b7af6..ed62767f 100644 --- a/src/geoptimiser.c +++ b/src/geoptimiser.c @@ -1,15 +1,15 @@ /* * geoptimiser.c * - * Refines detector geometry + * Refine detector geometry * - * Copyright © 2014 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2014 Thomas White <taw@physics.org> - * Oleksandr Yefanov - * Valerio Mariani + * 2014-2015 Oleksandr Yefanov + * 2014-2015 Valerio Mariani + * 2014-2015 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -27,7 +27,9 @@ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. * */ - +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif #include <stdlib.h> #include <stdio.h> @@ -38,18 +40,31 @@ #include <time.h> #include <float.h> +#ifdef HAVE_CAIRO +#ifdef HAVE_GTK +#define HAVE_SAVE_TO_PNG 1 +#include <cairo.h> +#include <gdk/gdk.h> +#endif /* HAVE_CAIRO */ +#endif /* HAVE_GTK */ + #include <detector.h> #include <stream.h> #include <version.h> #include <crystal.h> #include <image.h> #include <utils.h> +#include <render.h> + +#include "hdfsee-render.h" struct imagefeature; static void show_help(const char *s) { - printf("Syntax: %s [options] input.stream\n\n", s); + printf("Syntax: %s -i input.stream -g input.geom -o refined.geom " + "-c connected_rgcollection -q quadrant_rgcollection [options]\n", + s); printf( "Refines detector geometry.\n" "\n" @@ -64,6 +79,7 @@ static void show_help(const char *s) " -q, --quadrants=<rg_coll> Rigid group collection for quadrants.\n" " -c, --connected=<rg_coll> Rigid group collection for connected\n" " ASICs.\n" +" --no-error-maps Do not generate error map PNGs.\n" " -x, --min-num-peaks-per-pixel=<num> Minimum number of peaks per pixel.\n" " Default: 3. \n" " -p, --min-num-peaks-per-panel=<num> Minimum number of peaks per pixel.\n" @@ -75,8 +91,8 @@ static void show_help(const char *s) " --no-stretch Do not optimize distance offset.\n" " Default: distance offset is optimized\n" " -m --max-peak-dist=<num> Maximum distance between predicted and\n" -" detected peaks\n" -" Default: 4.0.\n" +" detected peaks (in pixels)\n" +" Default: 4.0 pixels.\n" ); } @@ -220,6 +236,28 @@ static double compute_average_clen (struct detector *det, char **clen_from, } +static double extract_f_from_stuff(const char *field_name, + struct stuff_from_stream* stuff) +{ + int i; + + char field_name_plus_equal[256]; + snprintf(field_name_plus_equal, 256, "hdf5%s = ", field_name); + + for ( i=0; i<stuff->n_fields; i++ ) { + + if ( strncmp(stuff->fields[i], field_name_plus_equal, + strlen(field_name_plus_equal)) == 0 ) { + return atoi(stuff->fields[i]+ + strlen(field_name_plus_equal)); + } + } + + ERROR("Failed to recover camera length from stream file\n"); + return -1; +} + + static struct pattern_list *read_patterns_from_steam_file(const char *infile, struct detector *det) { @@ -284,7 +322,8 @@ static struct pattern_list *read_patterns_from_steam_file(const char *infile, struct pattern **patterns_new; patterns_new = realloc(pattern_list->patterns, - (max_patterns+1024)*sizeof(struct pattern *)); + (max_patterns+1024)* + sizeof(struct pattern *)); if ( patterns_new == NULL ) { ERROR("Failed to allocate " "memory for loaded patterns.\n"); @@ -299,7 +338,8 @@ static struct pattern_list *read_patterns_from_steam_file(const char *infile, patn = malloc(sizeof(struct pattern)); if ( patn == NULL ) { - ERROR("Failed to allocate memory for loaded patterns.\n"); + ERROR("Failed to allocate memory for loaded " + "patterns.\n"); free(pattern_list->patterns); free(pattern_list); return NULL; @@ -314,7 +354,7 @@ static struct pattern_list *read_patterns_from_steam_file(const char *infile, avg_clen = compute_average_clen(det, &clen_from, &offset); if ( avg_clen == -1 ) { avg_clen = extract_f_from_stuff(clen_from, - cur.stuff_from_stream)*1e-3; + cur.stuff_from_stream)*1e-3; avg_clen += offset; } @@ -331,12 +371,14 @@ static struct pattern_list *read_patterns_from_steam_file(const char *infile, UnitCell *cell; UnitCell **new_unit_cells; - cell = crystal_get_cell(cur.crystals[0]); + cell = crystal_get_cell(cur.crystals[i]); new_unit_cells = realloc(patn->unit_cells, - (patn->n_unit_cells+1)*sizeof(UnitCell *)); + (patn->n_unit_cells+1)* + sizeof(UnitCell *)); if ( new_unit_cells == NULL ) { - ERROR("Failed to allocate memory for loaded patterns.\n"); + ERROR("Failed to allocate memory for " + "loaded patterns.\n"); free(pattern_list->patterns); free(pattern_list); free(patn); @@ -347,7 +389,8 @@ static struct pattern_list *read_patterns_from_steam_file(const char *infile, patn->n_unit_cells++; patn->unit_cells = new_unit_cells; - crystal_reflist = crystal_get_reflections(cur.crystals[i]); + crystal_reflist = crystal_get_reflections( + cur.crystals[i]); for ( refl = first_refl(crystal_reflist, &iter); refl != NULL; @@ -368,7 +411,8 @@ static struct pattern_list *read_patterns_from_steam_file(const char *infile, pattern_list->n_patterns++; if ( pattern_list->n_patterns%1000 == 0 ) { - STATUS("Loaded %i indexed patterns from %i total patterns\n", + STATUS("Loaded %i indexed patterns from %i " + "total patterns.\n", pattern_list->n_patterns, ++n_chunks); } @@ -378,7 +422,7 @@ static struct pattern_list *read_patterns_from_steam_file(const char *infile, close_stream(st); - STATUS("Found %d indexed patterns in file %s (from a total of %d)\n", + STATUS("Found %d indexed patterns in file %s (from a total of %d).\n", pattern_list->n_patterns, infile, n_chunks ); return pattern_list; @@ -424,12 +468,13 @@ static void compute_avg_cell_parameters(struct pattern_list *pattern_list, for ( cri=0; cri<ptn->n_unit_cells; cri++ ) { - cell_get_parameters(ptn->unit_cells[cri], &cpar[0], // a - &cpar[1], // b - &cpar[2], // c - &cpar[3], // alpha - &cpar[4], // beta - &cpar[5]); // gamma + cell_get_parameters(ptn->unit_cells[cri], + &cpar[0], // a + &cpar[1], // b + &cpar[2], // c + &cpar[3], // alpha + &cpar[4], // beta + &cpar[5]); // gamma cpar[0] *= 1e9; cpar[1] *= 1e9; @@ -454,16 +499,16 @@ static void compute_avg_cell_parameters(struct pattern_list *pattern_list, } STATUS("Average cell coordinates:\n"); - STATUS("Average a, b, c (in A): %6.3f, %6.3f, %6.3f\n", + STATUS("Average a, b, c (in nm): %6.3f, %6.3f, %6.3f\n", avcp[0],avcp[1],avcp[2]); STATUS("Minimum -Maximum a, b, c:\n" "\t%6.3f - %6.3f,\n" "\t%6.3f - %6.3f,\n" "\t%6.3f - %6.3f\n", minc[0], maxc[0], minc[1], maxc[1], minc[2], maxc[2]); - STATUS("Average alpha,beta,gamma: %6.3f, %6.3f, %6.3f\n", + STATUS("Average alpha,beta,gamma in degrees: %6.3f, %6.3f, %6.3f\n", avcp[3], avcp[4], avcp[5]); - STATUS("Minimum - Maximum alpha,beta,gamma:\n" + STATUS("Minimum - Maximum alpha,beta,gamma in degrees:\n" "\t%5.2f - %5.2f,\n" "\t%5.2f - %5.2f,\n" "\t%5.2f - %5.2f\n", @@ -518,7 +563,8 @@ static double compute_clen_to_use(struct pattern_list *pattern_list, int found = 0; for ( i=0; i<num_clens; i++ ) { - if ( fabs(pattern_list->patterns[cp]->clen-clens[i])<0.0001 ) { + if ( fabs(pattern_list->patterns[cp]->clen-clens[i]) + <0.0001 ) { clens_population[i]++; lambdas[i] += pattern_list->patterns[cp]->lambda; found = 1; @@ -535,14 +581,14 @@ static double compute_clen_to_use(struct pattern_list *pattern_list, double *lambdas_new; clens_population_new = realloc(clens_population, - (max_clens+1024)*sizeof(int)); + (max_clens+1024)*sizeof(int)); clens_new = realloc(clens_population, - (max_clens+1024)*sizeof(double)); + (max_clens+1024)*sizeof(double)); lambdas_new = realloc(clens_population, - (max_clens+1024)*sizeof(double)); + (max_clens+1024)*sizeof(double)); - if ( clens_new == NULL || clens_population_new == NULL || - lambdas_new == NULL) { + if ( clens_new == NULL || clens_population_new == NULL + || lambdas_new == NULL) { ERROR("Failed to allocate memory for " "camera length list\n"); free(clens); @@ -570,7 +616,8 @@ static double compute_clen_to_use(struct pattern_list *pattern_list, } if ( num_clens == 1 ) { - STATUS("All patterns have the same camera length: %f\n", clens[0]); + STATUS("All patterns have the same camera length: %f m.\n", + clens[0]); } else { STATUS("%i different camera lengths were found for the input " "patterns:\n", num_clens); @@ -584,17 +631,20 @@ static double compute_clen_to_use(struct pattern_list *pattern_list, irecistep = 1/cqu.u; - min_braggp_dist = fmin(fmin(irecistep/avcp[0],irecistep/avcp[1]), + min_braggp_dist = fmin(fmin(irecistep/avcp[0], + irecistep/avcp[1]), irecistep/avcp[2]); - STATUS("Camera length %0.4f was found %i times.\n" - "Minimum inter-bragg peak distance (based on average cell " - "parameters): %0.1f pixels\n",clens[i], clens_population[i], + STATUS("Camera length %0.4f m was found %i times.\n" + "Minimum inter-bragg peak distance (based on " + "average cell parameters): %0.1f pixels.\n", + clens[i], clens_population[i], min_braggp_dist); if ( min_braggp_dist<1.2*max_peak_distance ) { - STATUS("WARNING: The distance between Bragg peaks is too " - "small: " - "%0.1f < 1.2*%0.1f\n", min_braggp_dist, - max_peak_distance); + STATUS("WARNING: The distance between Bragg " + "peaks is too small: " + "%0.1f < 1.2*%0.1f pixels.\n", + min_braggp_dist, + max_peak_distance); } if ( clens_population[i] > clens_population[best_clen] ) { best_clen = i; @@ -604,7 +654,7 @@ static double compute_clen_to_use(struct pattern_list *pattern_list, } if ( only_best_distance ) { - STATUS("Only %i patterns with CLEN=%0.4f will be used.\n", + STATUS("Only %i patterns with CLEN=%0.4f m will be used.\n", clens_population[best_clen], clen_to_use); } @@ -639,7 +689,8 @@ static double comp_median(double *arr, unsigned int n) return arr[median] ; } - /* Find median of low, middle and high items; swap into position low */ + // Find median of low, middle and high items; swap into position + // low middle = (low + high) / 2; if ( arr[middle]>arr[high] ) { A = arr[middle]; @@ -657,12 +708,14 @@ static double comp_median(double *arr, unsigned int n) arr[low] = A; } - /* Swap low item (now in position middle) into position (low+1) */ + // Swap low item (now in position middle) into position + // (low+1) A = arr[middle]; arr[middle] = arr[low+1]; arr[low+1] = A; - /* Nibble from each end towards middle, swapping items when stuck */ + // Nibble from each end towards middle, swapping items when + // stuck ll = low + 1; hh = high; while (1) { @@ -738,36 +791,150 @@ static void free_all_curr_pix_displ(struct single_pix_displ *all_pix_displ, } +static int fill_pixel_statistics(int *num_pix_displ, + struct single_pix_displ** curr_pix_displ, + struct single_pix_displ* all_pix_displ, + struct connected_data *conn_data, + int ifs, int iss, int di, int aw, int dfv, + double *displ_x, + double *displ_y, double *displ_abs) +{ + double *cPxAfs; + double *cPxAss; + int cnu = 0; + + cPxAfs = calloc(num_pix_displ[ifs+aw*iss], sizeof(double)); + if ( cPxAfs == NULL ) { + ERROR("Failed to allocate memory for pixel statistics.\n"); + return 1; + } + cPxAss = calloc(num_pix_displ[ifs+aw*iss], sizeof(double)); + if ( cPxAss == NULL ) { + ERROR("Failed to allocate memory for pixel statistics.\n"); + free(cPxAfs); + return 1; + } + + curr_pix_displ[ifs+aw*iss] = &all_pix_displ[ifs+aw*iss]; + + while (1) { + if (curr_pix_displ[ifs+aw*iss]->dfs == dfv) break; + cPxAfs[cnu] = curr_pix_displ[ifs+aw*iss]->dfs; + cPxAss[cnu] = curr_pix_displ[ifs+aw*iss]->dss; + cnu++; + if ( curr_pix_displ[ifs+aw*iss]->ne == NULL ) break; + curr_pix_displ[ifs+aw*iss] = curr_pix_displ[ifs+aw*iss]->ne; + } + + if ( cnu < 1 ) return 0; + + displ_x[ifs+aw*iss] = comp_median(cPxAfs, cnu); + displ_y[ifs+aw*iss] = comp_median(cPxAss, cnu); + displ_abs[ifs+aw*iss] = modulus2d(displ_x[ifs+aw*iss], + displ_y[ifs+aw*iss]); + conn_data[di].n_peaks_in_conn++; + + free(cPxAfs); + free(cPxAss); + + return 0; +} + + + +static int compute_panel_statistics(struct rg_collection *connected, + int *num_pix_displ, + struct single_pix_displ** curr_pix_displ, + struct single_pix_displ* all_pix_displ, + struct connected_data *conn_data, + int di, int ip, int np, + double dfv, + int aw, double *displ_x, + double *displ_y, double *displ_abs) +{ + struct panel *p; + int ifs, iss; + + p = connected->rigid_groups[di]->panels[ip]; + + for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { + for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { + if ( num_pix_displ[ifs+aw*iss]>=np ) { + + int ret; + + ret = fill_pixel_statistics(num_pix_displ, + curr_pix_displ, + all_pix_displ, + conn_data, + ifs, iss, di, aw, + dfv, displ_x, + displ_y, displ_abs); + + if ( ret == -2 ) { + break; + } else if ( ret != 0 ) { + return ret; + } + + } else { + displ_x[ifs+aw*iss] = dfv; + displ_y[ifs+aw*iss] = dfv; + displ_abs[ifs+aw*iss] = dfv; + } + + } + } + return 0; +} + + + +static int allocate_next_element(struct single_pix_displ** curr_pix_displ, + int ipx) +{ + curr_pix_displ[ipx]->ne = malloc(sizeof(struct single_pix_displ)); + if ( curr_pix_displ[ipx]->ne == NULL ) { + ERROR("Failed to allocate memory for pixel statistics.\n"); + return 1; + } + + curr_pix_displ[ipx] = curr_pix_displ[ipx]->ne; + + return 0; +} + + static int compute_pixel_statistics(struct pattern_list *pattern_list, - struct detector *det, - struct rg_collection *connected, - struct rg_collection *quadrants, - int num_pix_in_slab, - int max_peak_distance, int array_width, - double default_fill_value, - double min_num_peaks_per_pixel, - double min_num_peaks_per_panel, - int only_best_distance, - double clen_to_use, - double *slab_to_x, double *slab_to_y, - struct connected_data *conn_data, - double *displ_x, - double *displ_y, double *displ_abs, - struct single_pix_displ* all_pix_displ, - struct single_pix_displ** curr_pix_displ, - int *num_pix_displ) + struct detector *det, + struct rg_collection *connected, + struct rg_collection *quadrants, + int num_pix_in_slab, + int max_peak_distance, int aw, + double dfv, + double min_num_peaks_per_pixel, + double min_num_peaks_per_panel, + int only_best_distance, + double clen_to_use, + double *slab_to_x, double *slab_to_y, + struct connected_data *conn_data, + double *displ_x, + double *displ_y, double *displ_abs, + struct single_pix_displ* all_pix_displ, + struct single_pix_displ** curr_pix_displ, + int *num_pix_displ) { int ipx, cp, ich, di, ip, np; for (di=0; di<connected->n_rigid_groups; di++) { conn_data[di].num_quad = find_quad_for_connected( - connected->rigid_groups[di], - quadrants); + connected->rigid_groups[di], + quadrants); conn_data[di].cang = 0.0; conn_data[di].cstr = 1.0; - conn_data[di].sh_x = default_fill_value; - conn_data[di].sh_y = default_fill_value; + conn_data[di].sh_x = dfv; + conn_data[di].sh_y = dfv; conn_data[di].num_peaks_per_pixel = 1; conn_data[di].name = connected->rigid_groups[di]->name; conn_data[di].n_peaks_in_conn = 0; @@ -775,8 +942,8 @@ static int compute_pixel_statistics(struct pattern_list *pattern_list, for ( ipx=0; ipx<num_pix_in_slab; ipx++ ) { - all_pix_displ[ipx].dfs = default_fill_value; - all_pix_displ[ipx].dss = default_fill_value; + all_pix_displ[ipx].dfs = dfv; + all_pix_displ[ipx].dss = dfv; all_pix_displ[ipx].ne = NULL; curr_pix_displ[ipx] = &all_pix_displ[ipx]; num_pix_displ[ipx] = 0; @@ -787,7 +954,8 @@ static int compute_pixel_statistics(struct pattern_list *pattern_list, ImageFeatureList *flist = pattern_list->patterns[cp]->im_list; if ( only_best_distance ) { - if ( fabs(pattern_list->patterns[cp]->clen-clen_to_use)>0.0001 ) { + if ( fabs(pattern_list->patterns[cp]->clen-clen_to_use) > + 0.0001 ) { continue; } } @@ -807,25 +975,25 @@ static int compute_pixel_statistics(struct pattern_list *pattern_list, struct imagefeature *imfe = image_get_feature(flist, ich); compute_x_y(det, imfe->fs, imfe->ss, &fx, &fy); - refl = find_closest_reflection(rlist, fx, fy, det, &min_dist); + refl = find_closest_reflection(rlist, fx, fy, det, + &min_dist); if ( refl == NULL ) continue; - if ( min_dist<max_peak_distance ) { + if ( min_dist < max_peak_distance ) { - int ipx = ((int)rint(imfe->fs) + array_width* + int ipx = ((int)rint(imfe->fs) + aw* (int)rint(imfe->ss)); if ( num_pix_displ[ipx]>0 ) { - curr_pix_displ[ipx]->ne = malloc(sizeof( - struct single_pix_displ)); - if ( curr_pix_displ[ipx]->ne == NULL ) { - ERROR("Failed to allocate memory for pixel " - "statistics.\n"); - return 1; - } - curr_pix_displ[ipx] = curr_pix_displ[ipx]->ne; + int ret; + + ret = allocate_next_element(curr_pix_displ, + ipx); + + if ( ret != 0) return ret; + } get_detector_pos(refl, &rfs, &rss); @@ -845,78 +1013,29 @@ static int compute_pixel_statistics(struct pattern_list *pattern_list, continue; } - for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) { + for (ip=0; ip<connected->rigid_groups[di]->n_panels; + ip++) { - struct panel *p; - int ifs, iss; + int ret; - p = connected->rigid_groups[di]->panels[ip]; + ret = compute_panel_statistics(connected, + num_pix_displ, + curr_pix_displ, + all_pix_displ, + conn_data, di, ip, + np, + dfv, + aw, + displ_x, + displ_y, + displ_abs); - for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { - for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { - if ( num_pix_displ[ifs+array_width*iss]>=np ) { - - double *cPxAfs; - double *cPxAss; - int cnu = 0; - - cPxAfs = calloc(num_pix_displ[ifs+array_width*iss], - sizeof(double)); - if ( cPxAfs == NULL ) { - ERROR("Failed to allocate memory for " - "pixel statistics.\n"); - return 1; - } - cPxAss = calloc(num_pix_displ[ifs+array_width*iss], - sizeof(double)); - if ( cPxAss == NULL ) { - ERROR("Failed to allocate memory for " - "pixel statistics.\n"); - free(cPxAfs); - return 1; - } - - curr_pix_displ[ifs+array_width*iss] = - &all_pix_displ[ifs+array_width*iss]; - - while (1) { - if (curr_pix_displ[ifs+array_width*iss]->dfs - == default_fill_value) break; - cPxAfs[cnu] = - curr_pix_displ[ifs+array_width*iss]->dfs; - cPxAss[cnu] = - curr_pix_displ[ifs+array_width*iss]->dss; - cnu++; - if ( curr_pix_displ[ifs+array_width*iss]->ne == - NULL ) break; - curr_pix_displ[ifs+array_width*iss] = - curr_pix_displ[ifs+array_width*iss]->ne; - - } - - if ( cnu<1 ) continue; - - displ_x[ifs+array_width*iss] = - comp_median(cPxAfs, cnu); - displ_y[ifs+array_width*iss] = - comp_median(cPxAss, cnu); - displ_abs[ifs+array_width*iss] = modulus2d( - displ_x[ifs+array_width*iss], - displ_y[ifs+array_width*iss]); - conn_data[di].n_peaks_in_conn++; - - free(cPxAfs); - free(cPxAss); - - } else { - displ_x[ifs+array_width*iss] = default_fill_value; - displ_y[ifs+array_width*iss] = default_fill_value; - displ_abs[ifs+array_width*iss] = default_fill_value; - } - } - } + if ( ret !=0 ) return ret; } - if ( conn_data[di].n_peaks_in_conn>=min_num_peaks_per_panel ) { + + + if ( conn_data[di].n_peaks_in_conn >= + min_num_peaks_per_panel ) { conn_data[di].num_peaks_per_pixel = np; } } @@ -927,7 +1046,7 @@ static int compute_pixel_statistics(struct pattern_list *pattern_list, static double compute_error(struct rg_collection *connected, - int array_width, + int aw, struct connected_data *conn_data, int *num_pix_displ, double *displ_abs) @@ -949,13 +1068,13 @@ static double compute_error(struct rg_collection *connected, for (ifs=p->min_fs; ifs<p->max_fs+1; ifs++) { for (iss=p->min_ss; iss<p->max_ss+1; iss++) { - if ( num_pix_displ[ifs+array_width*iss]>= + if ( num_pix_displ[ifs+aw*iss]>= conn_data[di].num_peaks_per_pixel ) { double cer; - cer = displ_abs[ifs+array_width*iss]* - displ_abs[ifs+array_width*iss]; + cer = displ_abs[ifs+aw*iss]* + displ_abs[ifs+aw*iss]; connected_error += cer; num_connected_error++; total_error += cer; @@ -965,12 +1084,14 @@ static double compute_error(struct rg_collection *connected, } } - if ( num_connected_error>0 ) { + if ( num_connected_error > 0 ) { + connected_error /= (double)num_connected_error; connected_error = sqrt(connected_error); - STATUS("Error for connected group %s (%d peaks): " - "<delta^2> = %0.4f\n", conn_data[di].name, + STATUS("Error for connected group %s: %d pixels with " + "more than %d peaks: <delta^2> = %0.4f pixels.\n", + conn_data[di].name, num_connected_error, conn_data[di].num_peaks_per_pixel, connected_error); } @@ -987,7 +1108,7 @@ static double compute_error(struct rg_collection *connected, } -static void fill_coordinate_matrices(struct detector *det, int array_width, +static void fill_coordinate_matrices(struct detector *det, int aw, double *slab_to_x, double *slab_to_y) { int pi; @@ -1006,8 +1127,8 @@ static void fill_coordinate_matrices(struct detector *det, int array_width, compute_x_y(det, ifs, iss, &xs, &ys); - slab_to_x[iss*array_width+ifs] = xs; - slab_to_y[iss*array_width+ifs] = ys; + slab_to_x[iss*aw+ifs] = xs; + slab_to_y[iss*aw+ifs] = ys; } } @@ -1070,14 +1191,19 @@ static int correct_empty_panels(struct rg_collection *quadrants, if ( conn_data[di].n_peaks_in_conn<min_num_peaks_per_panel ) { if (aver_num_ang[conn_data[di].num_quad]>0) { - conn_data[di].cang = aver_ang[conn_data[di].num_quad]; - conn_data[di].cstr = aver_str[conn_data[di].num_quad]; - STATUS("Connected group %s has not enough peaks (%i). Using" - " average angle: %0.4f\n", conn_data[di].name, - conn_data[di].n_peaks_in_conn, conn_data[di].cang); + conn_data[di].cang = + aver_ang[conn_data[di].num_quad]; + conn_data[di].cstr = + aver_str[conn_data[di].num_quad]; + STATUS("Connected group %s has not enough peaks " + "(%i). Using average angle: %0.4f degrees\n", + conn_data[di].name, + conn_data[di].n_peaks_in_conn, + conn_data[di].cang); } else { - STATUS("Connected group %s has not enough peaks (%i). Left " - "unchanged\n", conn_data[di].name, + STATUS("Connected group %s has not enough peaks " + "(%i). Left unchanged\n", + conn_data[di].name, conn_data[di].n_peaks_in_conn); } } @@ -1109,15 +1235,19 @@ static void correct_angle_and_stretch(struct rg_collection *connected, p = connected->rigid_groups[di]->panels[ip]; newx = - p->fsx*cos(conn_data[di].cang)-p->fsy*sin(conn_data[di].cang); + p->fsx*cos(conn_data[di].cang)- + p->fsy*sin(conn_data[di].cang); newy = - p->fsx*sin(conn_data[di].cang)+p->fsy*cos(conn_data[di].cang); + p->fsx*sin(conn_data[di].cang)+ + p->fsy*cos(conn_data[di].cang); p->fsx = newx; p->fsy = newy; newx = - p->ssx*cos(conn_data[di].cang)-p->ssy*sin(conn_data[di].cang); + p->ssx*cos(conn_data[di].cang)- + p->ssy*sin(conn_data[di].cang); newy = - p->ssx*sin(conn_data[di].cang)+p->ssy*cos(conn_data[di].cang); + p->ssx*sin(conn_data[di].cang)+ + p->ssy*cos(conn_data[di].cang); p->ssx = newx; p->ssy = newy; } @@ -1131,8 +1261,8 @@ static void correct_angle_and_stretch(struct rg_collection *connected, for (pi=0; pi<det->n_panels; pi++) { det->panels[pi].coffset -= use_clen*(1.0-stretch_coeff); } - STATUS("Using a single offset distance for the whole detector: %f\n", - det->panels[0].coffset); + STATUS("Using a single offset distance for the whole detector: " + "%f m.\n", det->panels[0].coffset); for ( di=0; di<connected->n_rigid_groups; di++ ) { conn_data[di].cstr = stretch_coeff; @@ -1140,9 +1270,10 @@ static void correct_angle_and_stretch(struct rg_collection *connected, } else { - STATUS("Using individual distances for rigid panels\n"); + STATUS("Using individual distances for rigid panels.\n"); for ( di=0; di<connected->n_rigid_groups; di++ ) { - for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) { + for ( ip=0; ip<connected->rigid_groups[di]->n_panels; + ip++ ) { struct panel *p; @@ -1176,17 +1307,55 @@ static void shift_panels(struct rg_collection *connected, } else { struct panel *p0; - double connected_panel_dist; + double delta_x, delta_y; p0 = connected->rigid_groups[di]->panels[0]; - connected_panel_dist = modulus2d( - p->cnx-p0->cnx/conn_data[di].cstr, - p->cny-p0->cny/conn_data[di].cstr - ); + delta_x = (p->cnx-p0->cnx/conn_data[di].cstr); + delta_y = (p->cny-p0->cny/conn_data[di].cstr); - p->cnx = p0->cnx + connected_panel_dist*p0->fsx; - p->cny = p0->cny + connected_panel_dist*p0->fsy; + p->cnx = p0->cnx + delta_x * cos(conn_data[di].cang) + - delta_y * sin(conn_data[di].cang); + p->cny = p0->cny + delta_x * sin(conn_data[di].cang) + + delta_y * cos(conn_data[di].cang); + } + } + } +} + + +static void recompute_panel(struct connected_data *conn_data, int di, int ip, + struct rg_collection *connected, + double *slab_to_x, + double *slab_to_y, + double *recomputed_slab_to_x, + double *recomputed_slab_to_y, + double *displ_x, double *displ_y, + double stretch_coeff, int aw, int *num_pix_displ) +{ + + double c_stretch; + struct panel *p; + int ifs, iss; + + c_stretch = conn_data[di].cstr; + + if ( fabs(c_stretch)<FLT_EPSILON ) c_stretch = + stretch_coeff; + + p = connected->rigid_groups[di]->panels[ip]; + + for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { + for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { + recomputed_slab_to_x[ifs+aw*iss] /= c_stretch; + recomputed_slab_to_y[ifs+aw*iss] /= c_stretch; + if ( num_pix_displ[ifs+aw*iss] >= + conn_data[di].num_peaks_per_pixel) { + + displ_x[ifs+aw*iss] -= (slab_to_x[ifs+aw*iss]- + recomputed_slab_to_x[ifs+aw*iss]); + displ_y[ifs+aw*iss] -= (slab_to_y[ifs+aw*iss]- + recomputed_slab_to_y[ifs+aw*iss]); } } } @@ -1198,7 +1367,7 @@ static void recompute_differences(struct rg_collection *connected, double *recomputed_slab_to_x, double *recomputed_slab_to_y, struct connected_data *conn_data, - int stretch_coeff, int array_width, + double stretch_coeff, int aw, double *displ_x, double *displ_y, int *num_pix_displ) { @@ -1208,45 +1377,76 @@ static void recompute_differences(struct rg_collection *connected, for ( di=0; di<connected->n_rigid_groups; di++ ) { for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) { - double c_stretch; - struct panel *p; - int ifs, iss; - c_stretch = conn_data[di].cstr; + recompute_panel(conn_data, di, ip, connected, + slab_to_x, slab_to_y, + recomputed_slab_to_x, + recomputed_slab_to_y, + displ_x, displ_y, + stretch_coeff, aw, num_pix_displ); + } + } +} - if ( fabs(c_stretch)<FLT_EPSILON ) c_stretch = stretch_coeff; - p = connected->rigid_groups[di]->panels[ip]; +static void fill_av_in_panel(struct rg_collection *connected, int di, int ip, + struct connected_data *conn_data, + int *num_pix_displ, int aw, + double *av_in_panel_fs, + double *av_in_panel_ss, + double *displ_x, double *displ_y) +{ + struct panel *p; + int ifs, iss; - for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { - for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { - recomputed_slab_to_x[ifs+array_width*iss] /= c_stretch; - recomputed_slab_to_y[ifs+array_width*iss] /= c_stretch; - if ( num_pix_displ[ifs+array_width*iss] >= - conn_data[di].num_peaks_per_pixel) { - - displ_x[ifs+array_width*iss] -= - (slab_to_x[ifs+array_width*iss]- - recomputed_slab_to_x[ifs+array_width*iss]); - displ_y[ifs+array_width*iss] -= - (slab_to_y[ifs+array_width*iss]- - recomputed_slab_to_y[ifs+array_width*iss]); - } - } + p = connected->rigid_groups[di]->panels[ip]; + + for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { + for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { + + if (num_pix_displ[ifs+aw*iss]>= + conn_data[di].num_peaks_per_pixel) { + av_in_panel_fs[conn_data[di].n_peaks_in_conn] = + displ_x[ifs+aw*iss]; + av_in_panel_ss[conn_data[di].n_peaks_in_conn] = + displ_y[ifs+aw*iss]; + conn_data[di].n_peaks_in_conn++; } } } } +static void fill_con_data_sh(struct connected_data *conn_data, + double *av_in_panel_fs, + double *av_in_panel_ss, int di, + double max_peak_distance) +{ + conn_data[di].sh_x = comp_median(av_in_panel_fs, + conn_data[di].n_peaks_in_conn); + conn_data[di].sh_y = comp_median(av_in_panel_ss, + conn_data[di].n_peaks_in_conn); + STATUS("Panel %s, num pixels: %i, shifts (in pixels) X,Y: %0.8f, %0.8f\n", + conn_data[di].name, conn_data[di].n_peaks_in_conn, + conn_data[di].sh_x, conn_data[di].sh_y); + if ( modulus2d(conn_data[di].sh_x, conn_data[di].sh_y) > + 0.8*max_peak_distance ) { + STATUS(" WARNING: absolute shift is: %0.1f > 0.8*%0.1f pixels." + " Increase the value of the max_peak_distance parameter!\n", + modulus2d(conn_data[di].sh_x, conn_data[di].sh_y), + max_peak_distance); + } +} + + static int compute_shifts(struct rg_collection *connected, struct connected_data *conn_data, - int *num_pix_displ, int array_width, + int *num_pix_displ, int aw, int min_num_peaks_per_panel, - double default_fill_value, double max_peak_distance, - double *displ_x, double *displ_y ) + double dfv, double max_peak_distance, + double *displ_x, double *displ_y) { - STATUS("Median for panels\n"); + STATUS("Median for panels.\n"); int di, ip; @@ -1281,45 +1481,21 @@ static int compute_shifts(struct rg_collection *connected, for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) { - struct panel *p; - int ifs, iss; + fill_av_in_panel(connected, di, ip, conn_data, + num_pix_displ, aw, av_in_panel_fs, + av_in_panel_ss, displ_x, displ_y); - p = connected->rigid_groups[di]->panels[ip]; - - for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { - for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { - - if (num_pix_displ[ifs+array_width*iss]>= - conn_data[di].num_peaks_per_pixel) { - av_in_panel_fs[conn_data[di].n_peaks_in_conn] = - displ_x[ifs+array_width*iss]; - av_in_panel_ss[conn_data[di].n_peaks_in_conn] = - displ_y[ifs+array_width*iss]; - conn_data[di].n_peaks_in_conn++; - } - } - } } if ( conn_data[di].n_peaks_in_conn>=min_num_peaks_per_panel ) { - conn_data[di].sh_x = - comp_median(av_in_panel_fs, conn_data[di].n_peaks_in_conn); - conn_data[di].sh_y = - comp_median(av_in_panel_ss, conn_data[di].n_peaks_in_conn); - STATUS("Panel %s, num pixels: %i, shifts X,Y: %0.8f, %0.8f\n", - conn_data[di].name, conn_data[di].n_peaks_in_conn, - conn_data[di].sh_x, conn_data[di].sh_y); - if ( modulus2d(conn_data[di].sh_x, conn_data[di].sh_y)> - 0.8*max_peak_distance ) { - STATUS(" WARNING: absolute shift is: %0.1f > 0.8*%0.1f. " - "Increase the value of the max_peak_distance parameter!\n", - modulus2d(conn_data[di].sh_x, conn_data[di].sh_y), - max_peak_distance); - } + fill_con_data_sh(conn_data, av_in_panel_fs, + av_in_panel_ss, di, + max_peak_distance); + } else { - conn_data[di].sh_x = default_fill_value; - conn_data[di].sh_y = default_fill_value; + conn_data[di].sh_x = dfv; + conn_data[di].sh_y = dfv; } free(av_in_panel_fs); free(av_in_panel_ss); @@ -1389,13 +1565,16 @@ static int compute_shifts_for_empty_panels(struct rg_collection *quadrants, if ( num_aver[conn_data[di].num_quad]>0 ) { conn_data[di].sh_x = aver_x[conn_data[di].num_quad]; conn_data[di].sh_y = aver_y[conn_data[di].num_quad]; - STATUS("Panel %s has not enough (%i) peaks. Using average " - "shifts X,Y: %0.2f,%0.2f\n", conn_data[di].name, + STATUS("Panel %s has not enough (%i) peaks. " + "Using average shifts (in pixels) X,Y: " + "%0.2f,%0.2f\n", conn_data[di].name, conn_data[di].n_peaks_in_conn, conn_data[di].sh_x, conn_data[di].sh_y); } else { - STATUS("Panel %s has not enough (%i) peaks. Left unchanged\n", - conn_data[di].name, conn_data[di].n_peaks_in_conn); + STATUS("Panel %s has not enough (%i) peaks. " + "Left unchanged.\n", + conn_data[di].name, + conn_data[di].n_peaks_in_conn); } } } @@ -1410,7 +1589,7 @@ static int compute_shifts_for_empty_panels(struct rg_collection *quadrants, static void correct_shifts(struct rg_collection *connected, struct connected_data *conn_data, - double default_fill_value, double clen_to_use) + double dfv, double clen_to_use) { int di; @@ -1423,8 +1602,8 @@ static void correct_shifts(struct rg_collection *connected, p = connected->rigid_groups[di]->panels[ip]; - if ( conn_data[di].sh_x>default_fill_value+1.0 && - conn_data[di].sh_y > default_fill_value+1.0 ) { + if ( conn_data[di].sh_x>dfv+1.0 && + conn_data[di].sh_y > dfv+1.0 ) { p->cnx -= conn_data[di].sh_x; p->cny -= conn_data[di].sh_y; @@ -1438,15 +1617,176 @@ static void correct_shifts(struct rg_collection *connected, } -static int compute_angles_and_stretch - (struct rg_collection *connected, +static void a_s_counting_loop(int *num_pix_displ, int ifs, int iss, + int di, struct connected_data *conn_data, + int aw, double *slab_to_x, + double *slab_to_y, struct panel *p0, + struct panel *p1, double *displ_x, + double *displ_y, double minrad, + int *num_ang) +{ + + double coX, coY, cdX, cdY; + + if ( num_pix_displ[ifs+aw*iss]>= + conn_data[di].num_peaks_per_pixel ) { + + int ifs1, iss1; + int max_fs1_tmp = p0->max_fs; + int max_ss1_tmp = p0->max_ss; + + coX = slab_to_x[ifs+aw*iss]; + coY = slab_to_y[ifs+aw*iss]; + cdX = coX - displ_x[ifs+aw*iss]; + cdY = coY - displ_y[ifs+aw*iss]; + + for (ifs1=ifs+1; ifs1<max_fs1_tmp+1; ifs1++) { + + if ( ifs1 == max_fs1_tmp ) { + max_fs1_tmp = p1->max_fs; + } + + for (iss1=iss+1; iss1<max_ss1_tmp+1; iss1++) { + + if ( iss1 == max_ss1_tmp ) { + max_ss1_tmp = p1->max_ss; + } + + if ( num_pix_displ[ifs1+aw*iss1]>= + conn_data[di].num_peaks_per_pixel ) { + + double dist; + double coX1, coY1, cdX1, cdY1; + double len1, len2; + + dist = modulus2d(ifs-ifs1,iss-iss1); + if ( dist < minrad ) continue; + coX1 = slab_to_x[ifs1+aw*iss1]; + coY1 = slab_to_y[ifs1+aw*iss1]; + cdX1 = + coX1 - displ_x[ifs1+aw*iss1]; + cdY1 = + coY1 - displ_y[ifs1+aw*iss1]; + + len1 = modulus2d(coX1-coX, coY1-coY); + len2 = modulus2d(cdX1-cdX, cdY1-cdY); + if ( len1<FLT_EPSILON || + len2<FLT_EPSILON ) { + continue; + } + + *num_ang = *num_ang+1; + } + } + } + } +} + + +static int a_s_processing_loop(int *num_pix_displ, int ifs, int iss, + int di, struct connected_data *conn_data, + int aw, double *slab_to_x, + double *slab_to_y, struct panel *p0, + struct panel *p1, double *displ_x, + double *displ_y, double minrad, + int max_num_ang, int *num_ang, + double *angles, double *stretches) +{ + double coX, coY, cdX, cdY; + + if ( num_pix_displ[ifs+aw*iss]>= + conn_data[di].num_peaks_per_pixel ) { + + int ifs1, iss1; + int max_fs1_tmp = p0->max_fs; + int max_ss1_tmp = p0->max_ss; + + if ( *num_ang>=max_num_ang ) return -2; + + coX = slab_to_x[ifs+aw*iss]; + coY = slab_to_y[ifs+aw*iss]; + cdX = coX - displ_x[ifs+aw*iss]; + cdY = coY - displ_y[ifs+aw*iss]; + + for (ifs1=ifs+1; ifs1<max_fs1_tmp+1; ifs1++) { + + if ( ifs1 == max_fs1_tmp ) { + max_fs1_tmp = p1->max_fs; + } + + for (iss1=iss+1; iss1<max_ss1_tmp+1; iss1++) { + + if ( iss1 == max_ss1_tmp ) { + max_ss1_tmp = p1->max_ss; + } + + if ( num_pix_displ[ifs1+aw*iss1]>= + conn_data[di].num_peaks_per_pixel ) { + + double dist; + double coX1, coY1, cdX1, cdY1; + double len1, len2; + double scalM; + double multlen; + + if ( *num_ang>=max_num_ang ) return -2; + dist = modulus2d(ifs-ifs1,iss-iss1); + if (dist<minrad) return 0; + coX1 = slab_to_x[ifs1+aw*iss1]; + coY1 = slab_to_y[ifs1+aw*iss1]; + cdX1 = + coX1 - displ_x[ifs1+aw*iss1]; + cdY1 = + coY1 - displ_y[ifs1+aw*iss1]; + + len1 = modulus2d(coX1-coX, coY1-coY); + len2 = modulus2d(cdX1-cdX, cdY1-cdY); + scalM = (coX1-coX)*(cdX1-cdX)+ + (coY1-coY)*(cdY1-cdY)- + FLT_EPSILON; + if ( len1<FLT_EPSILON || + len2<FLT_EPSILON ) { + return 0; + } + + multlen = len1*len2; + if ( fabs(scalM)>=multlen ) { + angles[*num_ang] = 0.0; + } else { + + angles[*num_ang] = 1.0; + + angles[*num_ang] = + acos(scalM/multlen); + + if ((coY1-coY)*(cdX1-cdX)- + (coX1-coX)*(cdY1-cdY) < 0) { + angles[*num_ang] *= -1.; + } + + } + + stretches[*num_ang] = len1/len2; + + *num_ang = *num_ang+1; + } + } + } + } + return 0; +} + + + + +static int compute_angles_and_stretch(struct rg_collection *connected, struct connected_data *conn_data, int *num_pix_displ, double *slab_to_x, double *slab_to_y, double *displ_x, double *displ_y, - int array_width, + int aw, int min_num_peaks_per_panel, double dist_coeff_ang_str, int num_peaks_per_pixel, @@ -1461,18 +1801,21 @@ static int compute_angles_and_stretch csaa = malloc(sizeof(struct connected_stretch_and_angles)); if ( csaa == NULL ) { - ERROR("Failed to allocate memory to compute angles and stretch.\n"); + ERROR("Failed to allocate memory to compute angles and " + "stretch.\n"); return 1; } csaa->stretch_coeff = malloc(connected->n_rigid_groups*sizeof(double)); if ( csaa->stretch_coeff == NULL ) { - ERROR("Failed to allocate memory to compute angles and stretch.\n"); + ERROR("Failed to allocate memory to compute angles and " + "stretch.\n"); free(csaa); return 1; } csaa->num_angles = malloc(connected->n_rigid_groups*sizeof(unsigned int)); if ( csaa->num_angles == NULL ) { - ERROR("Failed to allocate memory to compute angles and stretch.\n"); + ERROR("Failed to allocate memory to compute angles and " + "stretch.\n"); free(csaa->stretch_coeff); free(csaa); return 1; @@ -1481,7 +1824,9 @@ static int compute_angles_and_stretch csaa->num_coeff=0; for ( di=0; di<connected->n_rigid_groups; di++ ) { - if ( conn_data[di].n_peaks_in_conn<min_num_peaks_per_panel ) continue; + if ( conn_data[di].n_peaks_in_conn<min_num_peaks_per_panel ) { + continue; + } unsigned int max_num_ang = 0; @@ -1509,9 +1854,11 @@ static int compute_angles_and_stretch struct panel *p0 = connected->rigid_groups[di]->panels[ip0]; - for ( ip1=0; ip1<connected->rigid_groups[di]->n_panels; ip1++ ) { + for ( ip1=0; ip1<connected->rigid_groups[di]->n_panels; + ip1++ ) { - struct panel *p1 = connected->rigid_groups [di]->panels[ip1]; + struct panel *p1 = + connected->rigid_groups [di]->panels[ip1]; int ifs, iss; int min_fs_tmp = p0->min_fs; @@ -1526,67 +1873,21 @@ static int compute_angles_and_stretch max_fs_tmp = p1->max_fs; } - for (iss=min_ss_tmp; iss<max_ss_tmp+1; iss++) { + for (iss=min_ss_tmp; iss<max_ss_tmp+1; + iss++) { if ( iss == max_ss_tmp ) { min_ss_tmp = p1->min_ss; max_ss_tmp = p1->max_ss; } - double coX, coY, cdX, cdY; - - if ( num_pix_displ[ifs+array_width*iss]>= - conn_data[di].num_peaks_per_pixel ) { - - int ifs1, iss1; - int max_fs1_tmp = p0->max_fs; - int max_ss1_tmp = p0->max_ss; - - coX = slab_to_x[ifs+array_width*iss]; - coY = slab_to_y[ifs+array_width*iss]; - cdX = coX - displ_x[ifs+array_width*iss]; - cdY = coY - displ_y[ifs+array_width*iss]; - - for (ifs1=ifs+1; ifs1<max_fs1_tmp+1; ifs1++) { - - if ( ifs1 == max_fs1_tmp ) { - max_fs1_tmp = p1->max_fs; - } - - for (iss1=iss+1; iss1<max_ss1_tmp+1; iss1++) { + a_s_counting_loop(num_pix_displ, + ifs, iss, di, conn_data, + aw, slab_to_x, slab_to_y, + p0, p1, displ_x, + displ_y, minrad, + &num_ang); - if ( iss1 == max_ss1_tmp ) { - max_ss1_tmp = p1->max_ss; - } - - if ( num_pix_displ[ifs1+array_width*iss1]>= - conn_data[di].num_peaks_per_pixel ) { - - double dist; - double coX1, coY1, cdX1, cdY1; - double len1, len2; - - dist = modulus2d(ifs-ifs1,iss-iss1); - if ( dist < minrad ) continue; - coX1 = slab_to_x[ifs1+array_width*iss1]; - coY1 = slab_to_y[ifs1+array_width*iss1]; - cdX1 = - coX1 - displ_x[ifs1+array_width*iss1]; - cdY1 = - coY1 - displ_y[ifs1+array_width*iss1]; - - len1 = modulus2d(coX1-coX, coY1-coY); - len2 = modulus2d(cdX1-cdX, cdY1-cdY); - if ( len1<FLT_EPSILON || - len2<FLT_EPSILON ) { - continue; - } - - num_ang++; - } - } - } - } } } } @@ -1598,7 +1899,8 @@ static int compute_angles_and_stretch angles = malloc(max_num_ang*sizeof(double)); if ( angles == NULL ) { - ERROR("Error in allocating memory for angle optimization\n"); + ERROR("Error in allocating memory for angle " + "optimization\n"); free(csaa->stretch_coeff); free(csaa->num_angles); free(csaa); @@ -1606,7 +1908,8 @@ static int compute_angles_and_stretch } stretches = malloc(max_num_ang*sizeof(double)); if ( stretches == NULL ) { - ERROR("Error in allocating memory for stretch optimization\n"); + ERROR("Error in allocating memory for stretch " + "optimization\n"); free(angles); free(csaa->stretch_coeff); free(csaa->num_angles); @@ -1620,9 +1923,11 @@ static int compute_angles_and_stretch struct panel *p0 = connected->rigid_groups[di]->panels[ip0]; - for ( ip1=0; ip1<connected->rigid_groups[di]->n_panels; ip1++ ) { + for ( ip1=0; ip1<connected->rigid_groups[di]->n_panels; + ip1++ ) { - struct panel *p1 = connected->rigid_groups [di]->panels[ip1]; + struct panel *p1 = + connected->rigid_groups [di]->panels[ip1]; int ifs, iss; int min_fs_tmp = p0->min_fs; @@ -1637,94 +1942,29 @@ static int compute_angles_and_stretch max_fs_tmp = p1->max_fs; } - for (iss=min_ss_tmp; iss<max_ss_tmp+1; iss++) { + for (iss=min_ss_tmp; iss<max_ss_tmp+1; + iss++) { + + int ret; if ( iss == max_ss_tmp ) { min_ss_tmp = p1->min_ss; max_ss_tmp = p1->max_ss; } - double coX, coY, cdX, cdY; - - if ( num_pix_displ[ifs+array_width*iss]>= - conn_data[di].num_peaks_per_pixel ) { - - int ifs1, iss1; - int max_fs1_tmp = p0->max_fs; - int max_ss1_tmp = p0->max_ss; - - if ( num_ang>=max_num_ang ) break; - - coX = slab_to_x[ifs+array_width*iss]; - coY = slab_to_y[ifs+array_width*iss]; - cdX = coX - displ_x[ifs+array_width*iss]; - cdY = coY - displ_y[ifs+array_width*iss]; - - for (ifs1=ifs+1; ifs1<max_fs1_tmp+1; ifs1++) { - - if ( ifs1 == max_fs1_tmp ) { - max_fs1_tmp = p1->max_fs; - } - - for (iss1=iss+1; iss1<max_ss1_tmp+1; iss1++) { - - if ( iss1 == max_ss1_tmp ) { - max_ss1_tmp = p1->max_ss; - } - - if ( num_pix_displ[ifs1+array_width*iss1]>= - conn_data[di].num_peaks_per_pixel ) { - - double dist; - double coX1, coY1, cdX1, cdY1; - double len1, len2; - double scalM; - double multlen; - - if ( num_ang>=max_num_ang ) break; - dist = modulus2d(ifs-ifs1,iss-iss1); - if (dist<minrad) continue; - coX1 = slab_to_x[ifs1+array_width*iss1]; - coY1 = slab_to_y[ifs1+array_width*iss1]; - cdX1 = - coX1 - displ_x[ifs1+array_width*iss1]; - cdY1 = - coY1 - displ_y[ifs1+array_width*iss1]; - - len1 = modulus2d(coX1-coX, coY1-coY); - len2 = modulus2d(cdX1-cdX, cdY1-cdY); - scalM = (coX1-coX)*(cdX1-cdX)+ - (coY1-coY)*(cdY1-cdY)- - FLT_EPSILON; - if ( len1<FLT_EPSILON || - len2<FLT_EPSILON ) { - continue; - } - - multlen = len1*len2; - if ( fabs(scalM)>=multlen ) { - angles[num_ang] = 0.0; - } else { - - angles[num_ang] = 1.0; - - angles[num_ang] = - acos(scalM/multlen); - - if ((coY1-coY)*(cdX1-cdX)- - (coX1-coX)*(cdY1-cdY) < 0) { - angles[num_ang] *= -1.; - } - - } - - stretches[num_ang] = len1/len2; - - num_ang++; - } - } - } - } + ret = a_s_processing_loop( + num_pix_displ, + ifs, iss, di, + conn_data, + aw, slab_to_x, + slab_to_y, + p0, p1, displ_x, + displ_y, minrad, + max_num_ang, + &num_ang, angles, + stretches); + + if ( ret == -2 ) break; } } } @@ -1734,7 +1974,8 @@ static int compute_angles_and_stretch conn_data[di].cang = -comp_median(angles,num_ang); conn_data[di].cstr = comp_median(stretches,num_ang); - STATUS("Panel %s, num: %i, angle: %0.4f, stretch: %0.4f\n", + STATUS("Panel %s, num: %i, angle: %0.4f deg, stretch coeff: " + "%0.4f\n", conn_data[di].name, num_ang, conn_data[di].cang, conn_data[di].cstr); @@ -1773,8 +2014,9 @@ static int compute_angles_and_stretch stretch_cf = 0; for ( di=0; di<num_coeff; di++ ) { if ( conn_data[di].num_peaks_per_pixel>=ipp ) { - stretch_cf += total_num*csaa->stretch_coeff[di]* - (double)csaa->num_angles[di]; + stretch_cf += + total_num*csaa->stretch_coeff[di]* + (double)csaa->num_angles[di]; } } break; @@ -1789,7 +2031,8 @@ static int compute_angles_and_stretch stretch_cf); if ( man_stretching_coeff>FLT_EPSILON ) { stretch_cf = man_stretching_coeff; - STATUS("Using manually set stretch coefficient: %0.4f\n", stretch_cf); + STATUS("Using manually set stretch coefficient: %0.4f\n", + stretch_cf); for ( di=0; di<connected->n_rigid_groups; di++ ) { conn_data[di].cstr = man_stretching_coeff; @@ -1807,83 +2050,371 @@ static int compute_angles_and_stretch } -static int save_data_to_hdf5(char * filename, struct detector* det, - int max_fs, int max_ss, double default_fill_value, - double *data) +#ifdef HAVE_SAVE_TO_PNG + +static void draw_panel(struct image *image, cairo_t *cr, + cairo_matrix_t *basic_m, GdkPixbuf **pixbufs, int i) { + struct panel p = image->det->panels[i]; + int w = gdk_pixbuf_get_width(pixbufs[i]); + int h = gdk_pixbuf_get_height(pixbufs[i]); + cairo_matrix_t m; + + /* Start with the basic coordinate system */ + cairo_set_matrix(cr, basic_m); + + /* Move to the right location */ + cairo_translate(cr, p.cnx, p.cny); + + /* Twiddle directions according to matrix */ + cairo_matrix_init(&m, p.fsx, p.fsy, p.ssx, p.ssy, 0.0, 0.0); + cairo_transform(cr, &m); + + gdk_cairo_set_source_pixbuf(cr, pixbufs[i], 0.0, 0.0); + cairo_rectangle(cr, 0.0, 0.0, w, h); +} + + +struct rectangle { - struct image *im; - int i; - int ret; + int width, height; + double min_x, min_y, max_x, max_y; +}; - im = malloc(sizeof(struct image)); - if ( im == NULL ) { - ERROR("Failed to allocate memory to save data.\n"); + +static int unpack_slab(struct image *image) +{ + struct detector *det = image->det; + int pi; + + image->dp = malloc(det->n_panels * sizeof(float *)); + image->bad = malloc(det->n_panels * sizeof(int *)); + if ( (image->dp == NULL) || (image->bad == NULL) ) { + ERROR("Failed to allocate panels.\n"); return 1; } - im->data = malloc((max_fs+1)*(max_ss+1)*sizeof(float)); - if ( im->data == NULL ) { + + for ( pi=0; pi<det->n_panels; pi++ ) { + + struct panel *p; + int fs, ss; + + p = &det->panels[pi]; + image->dp[pi] = malloc(p->w*p->h*sizeof(float)); + image->bad[pi] = calloc(p->w*p->h, sizeof(int)); + if ( (image->dp[pi] == NULL) || (image->bad[pi] == NULL) ) { + ERROR("Failed to allocate panel\n"); + return 1; + } + + for ( ss=0; ss<p->h; ss++ ) { + for ( fs=0; fs<p->w; fs++ ) { + + int idx; + int cfs, css; + + cfs = fs+p->min_fs; + css = ss+p->min_ss; + idx = cfs + css*image->width; + + image->dp[pi][fs+p->w*ss] = image->data[idx]; + image->bad[pi][fs+p->w*ss] = 0; + + } + } + } + + return 0; +} + + +static int draw_detector(cairo_surface_t *surf, struct image *image, + struct rectangle rect) +{ + cairo_t *cr; + cairo_matrix_t basic_m; + cairo_matrix_t m; + GdkPixbuf **pixbufs; + int n_pixbufs; + + cr = cairo_create(surf); + + unpack_slab(image); + pixbufs = render_panels(image, 1, SCALE_GEOPTIMISER, 1, &n_pixbufs); + + /* Blank grey background */ + cairo_rectangle(cr, 0.0, 0.0, rect.width, rect.height); + cairo_set_source_rgb(cr, 0.5, 0.5, 0.5); + cairo_fill(cr); + + /* Set up basic coordinate system + * - origin in the centre, y upwards. */ + cairo_identity_matrix(cr); + cairo_matrix_init(&m, 1.0, 0.0, 0.0, -1.0, 0.0, 0.0); + + + cairo_translate(cr, -rect.min_x , rect.max_y); + cairo_transform(cr, &m); + cairo_get_matrix(cr, &basic_m); + + if (pixbufs != NULL) { + + int i; + + for (i = 0; i < image->det->n_panels; i++) { + draw_panel(image, cr, &basic_m, pixbufs, i); + cairo_fill(cr); + } + + } + + /* Free old pixbufs */ + if (pixbufs != NULL) { + int i; + for (i = 0; i < n_pixbufs; i++) { + g_object_unref(pixbufs[i]); + } + free(pixbufs); + } + + return 0; + +} + + +static int save_data_to_png(char *filename, struct detector *det, + int max_fs, int max_ss, double default_fill_value, + double *data) +{ + struct image im; + int i; + struct rectangle rect; + GdkPixbuf *col_scale; + cairo_t *cr; + + cairo_status_t r; + cairo_surface_t *surf; + + im.data = malloc((max_fs+1)*(max_ss+1)*sizeof(float)); + if ( im.data == NULL ) { ERROR("Failed to allocate memory to save data.\n"); - free(im); return 1; } - im->det = det; - im->width = max_fs+1; - im->height = max_ss+1; - im->beam = NULL; - im->spectrum = NULL; + im.det = det; + im.width = max_fs+1; + im.height = max_ss+1; + im.flags = NULL; for ( i=0; i<(max_fs+1)*(max_ss+1); i++) { if ( data[i] == default_fill_value ) { - im->data[i] = 0.0; + im.data[i] = 0.0; + } else if ( data[i] > 1.0) { + im.data[i] = 1.0; } else { - im->data[i] = (float)data[i]; + im.data[i] = (float)data[i]; } + im.data[i] *= 10.0; /* render_panels sets this as max */ } - ret = hdf5_write_image(filename, im, NULL); + get_pixel_extents(im.det, &rect.min_x, &rect.min_y, &rect.max_x, + &rect.max_y); - if ( ret != 0 ) { - free(im->data); - free(im); + if (rect.min_x > 0.0) rect.min_x = 0.0; + if (rect.max_x < 0.0) rect.max_x = 0.0; + if (rect.min_y > 0.0) rect.min_y = 0.0; + if (rect.max_y < 0.0) rect.max_y = 0.0; + + rect.width = rect.max_x - rect.min_x; + rect.height = rect.max_y - rect.min_y; + + /* Add a thin border */ + rect.width += 2.0; + rect.height += 2.0; + surf = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, rect.width + 20, + rect.height); + + draw_detector(surf, &im, rect); + + col_scale = render_get_colour_scale(20, rect.height, SCALE_GEOPTIMISER); + + cr = cairo_create(surf); + cairo_identity_matrix(cr); + cairo_translate(cr, rect.width, 0.0); + cairo_rectangle(cr, 0.0, 0.0, 20.0, rect.height); + gdk_cairo_set_source_pixbuf(cr, col_scale, 0.0, 0.0); + cairo_fill(cr); + cairo_destroy(cr); + + r = cairo_surface_write_to_png(surf, filename); + if (r != CAIRO_STATUS_SUCCESS) { + free(im.data); return 1; } - free(im->data); - free(im); + free(im.data); return 0; } +#endif /* HAVE_SAVE_TO_PNG */ + +static void calculate_panel_correction(int di, int ip, int aw, + int *num_pix_displ, + struct rg_collection *connected, + struct connected_data *conn_data) +{ + struct panel *p; + int ifs, iss; + + p = connected->rigid_groups[di]->panels[ip]; + for (ifs=p->min_fs; ifs<p->max_fs+1; ifs++) { + for (iss=p->min_ss; iss<p->max_ss+1; iss++) { + if ( num_pix_displ[ifs+aw*iss]>= + conn_data[di].num_peaks_per_pixel ) { + conn_data[di].n_peaks_in_conn++; + } + } + } + +} + + +static void compute_abs_displ(struct rg_collection *connected, + struct connected_data *conn_data, + int *num_pix_displ, + double dfv, int di, int ip, int aw, + double *displ_x, + double *displ_y, + double *displ_abs) +{ + struct panel *p; + int ifs, iss; + + if (conn_data[di].sh_x < dfv+1) return; + + p = connected->rigid_groups[di]->panels[ip]; + + for (ifs=p->min_fs; ifs<p->max_fs+1; ifs++) { + for (iss=p->min_ss; iss<p->max_ss+1; iss++) { + if ( num_pix_displ[ifs+aw*iss]>= + conn_data[di].num_peaks_per_pixel ) { + displ_x[ifs+aw*iss] -= conn_data[di].sh_x; + displ_y[ifs+aw*iss] -= conn_data[di].sh_y; + displ_abs[ifs+aw*iss] = modulus2d( + displ_x[ifs+aw*iss], + displ_y[ifs+aw*iss] + ); + } else { + displ_abs[ifs+aw*iss] = dfv; + } + } + } +} + + +int check_and_enforce_cspad_dist(struct detector *det, int enforce) +{ + int np = 0; + int num_errors_found = 0; + + double dist_to_check = 197.0; + double tol = 0.2; + + for ( np=0; np<det->n_panels; np = np+2 ) { + + double dist2; + + struct panel *ep = &det->panels[np]; + struct panel *op = &det->panels[np+1]; + + dist2 = (( ep->cnx - op->cnx )*( ep->cnx - op->cnx ) + + ( ep->cny - op->cny )*( ep->cny - op->cny )); + + if ( dist2 > (dist_to_check+tol)*(dist_to_check+tol) || + dist2 < (dist_to_check-tol)*(dist_to_check-tol) ) { + + num_errors_found += 1; + + STATUS("Warning: distance between panels %s and %s " + "is outside acceptable margins.\n", ep->name, + op->name); + + if ( enforce ) { + + double new_op_cx, new_op_cy; + + new_op_cx = ep->cnx + ep->fsx*dist_to_check; + new_op_cy = ep->cny + ep->fsy*dist_to_check; + + op->cnx = new_op_cx; + op->cny = new_op_cy; + + STATUS("Enforcing distance....\n"); + } + + } + + if ( ep->fsx != op->fsx || ep->ssx != op->ssx || + ep->fsy != op->fsy || ep->ssx != op->ssx ) { + + num_errors_found += 1; + + STATUS("Warning: relative orientation between panels " + "%s and %s is incorrect.\n", ep->name, op->name); + + if ( enforce ) { + + STATUS("Enforcing relative orientation....\n"); + + op->fsx = ep->fsx; + op->ssx = ep->ssx; + op->fsy = ep->fsy; + op->ssy = ep->ssy; + + op->xfs = ep->xfs; + op->xss = ep->xss; + op->yfs = ep->yfs; + op->yss = ep->yss; + } + + } + + } + return num_errors_found; +} + + int optimize_geometry(char *infile, char *outfile, char *geometry_filename, struct detector *det, struct rg_collection* quadrants, struct rg_collection* connected, int min_num_peaks_per_pixel, int min_num_peaks_per_panel, int only_best_distance, int nostretch, - int individual_coffset, double max_peak_dist, - const char *command_line) + int individual_coffset, int error_maps, + int enforce_cspad_layout, int no_cspad, + double max_peak_dist, const char *command_line) { int num_pix_in_slab; int max_fs = 0; int max_ss = 0; - int array_width = 0; + int aw = 0; int pi, di, ip, pti; - int ret1, ret2, ret3; - int ret4, ret5, ret6; + int ret1, ret2; int ret; int write_ret; + int maybe_cspad = 0; + + int *num_pix_displ; double res_sum; double istep; double clen_to_use; double man_stretching_coeff = 0.0; double avc[6] = {0.,0.,0.,0.,0.,0.}; - double default_fill_value = -10000.0; - double dist_coeff_ang_str = 0.2; // for angles and stretch calculation use - // only pixels which are distco*size_panel - // away - - int *num_pix_displ; + double dfv = -10000.0; + // for angles and stretch calculation use + // only pixels which are distco*size_panel + // away + double dist_coeff_ang_str = 0.2; double *displ_x; double *displ_y; double *displ_abs; @@ -1894,6 +2425,7 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, double* recomputed_slab_to_x; double* recomputed_slab_to_y; double stretch_coeff = 1; + struct single_pix_displ *all_pix_displ; struct single_pix_displ **curr_pix_displ; struct connected_data *conn_data = NULL; @@ -1901,13 +2433,62 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, if ( nostretch ) man_stretching_coeff = 1.0; - STATUS("Maximum distance between peaks: %0.1f\n", max_peak_dist); + STATUS("Maximum distance between peaks: %0.1f pixels.\n", max_peak_dist); STATUS("Minimum number of measurements for pixel to be included in the " "refinement: %i\n", min_num_peaks_per_pixel); - STATUS("Minimum number of measurements for panel for accurate estimation of" - " position/orientation: %i\n", min_num_peaks_per_panel); + STATUS("Minimum number of measurements for panel for accurate estimation " + "of position/orientation: %i\n", min_num_peaks_per_panel); + + if ( det->n_panels == 64 ) { + maybe_cspad = 1; + } + + if ( maybe_cspad && !no_cspad ) { + + int num_errors = 0; + + STATUS("It looks like the detector is a CSPAD. " + "Checking relative distance and orientation of " + "connected ASICS.\n"); + STATUS("If the detector is not a CSPAD, please rerun the " + "program with the --no-cspad option.\n"); + if ( enforce_cspad_layout ) { + STATUS("Enforcing CSPAD layout...\n"); + } + + num_errors = check_and_enforce_cspad_dist(det, + enforce_cspad_layout); + + if ( enforce_cspad_layout ) { + + int geom_wr; + + STATUS("Saving geometry with enforced CSPAD layout.\n" + "Please restart geometry optimization using the " + "optimized geometry from this run as input geometry " + "file.\n"); + geom_wr = write_detector_geometry_2(geometry_filename, + outfile, det, + command_line, 1); + if ( geom_wr != 0 ) { + ERROR("Error in writing output geometry file.\n"); + return 1; + } + STATUS("All done!\n"); + return 0; + } + + if ( !enforce_cspad_layout && num_errors > 0 ) { + ERROR("Relative distance and orientation of connected " + "ASICS do not respect the CSPAD layout.\n" + "Geometry optimization cannot continue.\n" + "Please rerun the program with the " + "--enforce-cspad-layout option.\n"); + return 1; + } + } pattern_list = read_patterns_from_steam_file(infile, det); if ( pattern_list->n_patterns < 1 ) { @@ -1931,7 +2512,7 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, istep = res_sum/det->n_panels; - array_width = max_fs+1; + aw = max_fs+1; clen_to_use = compute_clen_to_use(pattern_list, istep, avc, max_peak_dist, @@ -1978,9 +2559,10 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, return 1; } - fill_coordinate_matrices(det, array_width, slab_to_x, slab_to_y); + fill_coordinate_matrices(det, aw, slab_to_x, slab_to_y); - all_pix_displ = calloc(num_pix_in_slab, sizeof(struct single_pix_displ)); + all_pix_displ = calloc(num_pix_in_slab, + sizeof(struct single_pix_displ)); if ( all_pix_displ == NULL ) { ERROR("Error allocating memory for connected structure data.\n"); free(displ_x); @@ -1990,7 +2572,8 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, free(slab_to_y); return 1; } - curr_pix_displ = calloc(num_pix_in_slab, sizeof(struct single_pix_displ*)); + curr_pix_displ = calloc(num_pix_in_slab, + sizeof(struct single_pix_displ*)); if ( curr_pix_displ == NULL ) { ERROR("Error allocating memory for connected structure data.\n"); free(displ_x); @@ -2014,7 +2597,8 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, return 1; } - conn_data = malloc(connected->n_rigid_groups*sizeof(struct connected_data)); + conn_data = malloc(connected->n_rigid_groups* + sizeof(struct connected_data)); if ( conn_data == NULL ) { ERROR("Error allocating memory for connected structure data.\n"); free(displ_x); @@ -2028,14 +2612,17 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, return 1; } - STATUS("Computing pixel statistics\n"); + STATUS("Computing pixel statistics.\n"); ret = compute_pixel_statistics(pattern_list, det, connected, quadrants, num_pix_in_slab, max_peak_dist, - array_width, default_fill_value, + aw, dfv, min_num_peaks_per_pixel, - min_num_peaks_per_panel, only_best_distance, - clen_to_use, slab_to_x, slab_to_y, conn_data, - displ_x, displ_y, displ_abs, all_pix_displ, + min_num_peaks_per_panel, + only_best_distance, + clen_to_use, slab_to_x, + slab_to_y, conn_data, + displ_x, displ_y, displ_abs, + all_pix_displ, curr_pix_displ, num_pix_displ); if ( ret != 0 ) { @@ -2044,7 +2631,8 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, free(displ_abs); free(slab_to_x); free(slab_to_y); - free_all_curr_pix_displ(all_pix_displ, curr_pix_displ, num_pix_in_slab); + free_all_curr_pix_displ(all_pix_displ, curr_pix_displ, + num_pix_in_slab); free(num_pix_displ); free(conn_data); return 1; @@ -2056,7 +2644,8 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, image_feature_list_free(pattern_list->patterns[pti]->im_list); reflist_free(pattern_list->patterns[pti]->ref_list); - for ( nuc=0; nuc<pattern_list->patterns[pti]->n_unit_cells; nuc++) { + for ( nuc=0; nuc<pattern_list->patterns[pti]->n_unit_cells; + nuc++) { cell_free(pattern_list->patterns[pti]->unit_cells[nuc]); } free(pattern_list->patterns[pti]->filename); @@ -2064,31 +2653,41 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, } free(pattern_list); - STATUS("Saving displacements before corrections\n"); - ret1 = save_data_to_hdf5("disp_x_before.h5", det, max_fs, max_ss, - default_fill_value, displ_x); - ret2 = save_data_to_hdf5("disp_y_before.h5", det, max_fs, max_ss, - default_fill_value, displ_y); - ret3 = save_data_to_hdf5("disp_abs_before.h5", det, max_fs, max_ss, - default_fill_value, displ_abs); - if ( ret1!=0 || ret2!=0 || ret3!=0 ) { - ERROR("Error while writing data to file.\n"); - free(conn_data); - free(displ_x); - free(displ_y); - free(displ_abs); - free(num_pix_displ); - free(slab_to_x); - free(slab_to_y); - return 1; + if ( error_maps ) { + STATUS("Saving displacements before corrections.\n"); + +#ifdef HAVE_SAVE_TO_PNG + + ret1 = save_data_to_png("error_map_before.png", det, max_fs, max_ss, + dfv, displ_abs); + if ( ret1!=0 ) { + ERROR("Error while writing data to file.\n"); + free(conn_data); + free(displ_x); + free(displ_y); + free(displ_abs); + free(num_pix_displ); + free(slab_to_x); + free(slab_to_y); + return 1; + } + +#else /* HAVE_SAVE_TO_PNG */ + + STATUS("ERROR: geoptimiser was compiled without GTK and cairo " + "support. Error maps will not be saved.\n"); + +#endif /* HAVE_SAVE_TO_PNG */ + } STATUS("Computing initial error.\n"); - totalError = compute_error(connected, array_width, conn_data, + totalError = compute_error(connected, aw, conn_data, num_pix_displ, displ_abs); - STATUS("The total initial error <delta^2> = %0.4f\n", totalError); - STATUS("Now calculating corrections\n"); + STATUS("The total initial error <delta^2> = %0.4f pixels.\n", + totalError); + STATUS("Now calculating corrections.\n"); for ( di=0;di<connected->n_rigid_groups;di++ ) { @@ -2096,28 +2695,19 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) { - struct panel *p; - int ifs, iss; + calculate_panel_correction(di, ip, aw, num_pix_displ, + connected, conn_data); - p = connected->rigid_groups[di]->panels[ip]; - for (ifs=p->min_fs; ifs<p->max_fs+1; ifs++) { - for (iss=p->min_ss; iss<p->max_ss+1; iss++) { - if ( num_pix_displ[ifs+array_width*iss]>= - conn_data[di].num_peaks_per_pixel ) { - conn_data[di].n_peaks_in_conn++; - } - } - } } } - STATUS("Calculating angles and elongations (usually long)\n"); + STATUS("Calculating angles and elongations.\n"); ret = compute_angles_and_stretch(connected, conn_data, num_pix_displ, slab_to_x, slab_to_y, displ_x, displ_y, - array_width, + aw, min_num_peaks_per_panel, dist_coeff_ang_str, min_num_peaks_per_pixel, @@ -2178,16 +2768,17 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, return 1; } - fill_coordinate_matrices(det, array_width, recomputed_slab_to_x, + fill_coordinate_matrices(det, aw, recomputed_slab_to_x, recomputed_slab_to_y); - recompute_differences(connected, slab_to_x, slab_to_y, recomputed_slab_to_x, + recompute_differences(connected, slab_to_x, slab_to_y, + recomputed_slab_to_x, recomputed_slab_to_y, conn_data, - stretch_coeff, array_width, displ_x, displ_y, + stretch_coeff, aw, displ_x, displ_y, num_pix_displ); - ret = compute_shifts(connected, conn_data, num_pix_displ, array_width, - min_num_peaks_per_panel, default_fill_value, + ret = compute_shifts(connected, conn_data, num_pix_displ, aw, + min_num_peaks_per_panel, dfv, max_peak_dist, displ_x, displ_y ); if ( ret != 0 ) return 1; @@ -2198,59 +2789,51 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, for ( di=0;di<connected->n_rigid_groups;di++ ) { for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) { - struct panel *p; - int ifs, iss; + compute_abs_displ(connected, conn_data, + num_pix_displ, dfv, di, ip, aw, + displ_x, displ_y, displ_abs); - if (conn_data[di].sh_x < default_fill_value+1) continue; + } + } - p = connected->rigid_groups[di]->panels[ip]; + correct_shifts(connected, conn_data, dfv, clen_to_use); - for (ifs=p->min_fs; ifs<p->max_fs+1; ifs++) { - for (iss=p->min_ss; iss<p->max_ss+1; iss++) { - if ( num_pix_displ[ifs+array_width*iss]>= - conn_data[di].num_peaks_per_pixel ) { - displ_x[ifs+array_width*iss] -= conn_data[di].sh_x; - displ_y[ifs+array_width*iss] -= conn_data[di].sh_y; - displ_abs[ifs+array_width*iss] = modulus2d( - displ_x[ifs+array_width*iss], - displ_y[ifs+array_width*iss] - ); - } else { - displ_abs[ifs+array_width*iss] = default_fill_value; - } - } - } + if ( error_maps ) { + + +#ifdef HAVE_SAVE_TO_PNG + + STATUS("Saving displacements after corrections.\n"); + ret2 = save_data_to_png("error_map_after.png", det, max_fs, max_ss, + dfv, displ_x); + if ( ret2 !=0 ) { + ERROR("Error while writing data to file.\n"); + free(conn_data); + free(displ_x); + free(displ_y); + free(displ_abs); + free(num_pix_displ); + free(slab_to_x); + free(slab_to_y); + free(recomputed_slab_to_x); + free(recomputed_slab_to_y); + return 1; } - } - correct_shifts(connected, conn_data, default_fill_value, clen_to_use); +#else /* HAVE_SAVE_TO_PNG */ + + STATUS("ERROR: geoptimiser was compiled without GTK and cairo support.\n" + "Error maps will not be saved.\n"); + +#endif /* HAVE_SAVE_TO_PNG */ - STATUS("Saving displacements after corrections\n"); - ret4 = save_data_to_hdf5("disp_x_after.h5", det, max_fs, max_ss, - default_fill_value, displ_x); - ret5 = save_data_to_hdf5("disp_y_after.h5", det, max_fs, max_ss, - default_fill_value, displ_y); - ret6 = save_data_to_hdf5("disp_abs_after.h5", det, max_fs, max_ss, - default_fill_value, displ_abs); - if ( ret4!=0 || ret5!=0 || ret6!=0 ) { - ERROR("Error while writing data to file.\n"); - free(conn_data); - free(displ_x); - free(displ_y); - free(displ_abs); - free(num_pix_displ); - free(slab_to_x); - free(slab_to_y); - free(recomputed_slab_to_x); - free(recomputed_slab_to_y); - return 1; } STATUS("Computing final error.\n"); - totalError = compute_error(connected, array_width, conn_data, num_pix_displ, + totalError = compute_error(connected, aw, conn_data, num_pix_displ, displ_abs); - STATUS("The total final error <delta^2> = %0.4f\n",totalError); + STATUS("The total final error <delta^2> = %0.4f pixels.\n",totalError); write_ret = write_detector_geometry_2(geometry_filename, outfile, det, command_line, 1); @@ -2259,6 +2842,15 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, return 1; } STATUS("All done!\n"); + if ( error_maps ) { + +#ifdef HAVE_SAVE_TO_PNG + + STATUS("Be sure to inspect error_map_before.png and " + "error_map_after.png !!\n"); + +#endif /* HAVE_SAVE_TO_PNG */ + } free(conn_data); free(displ_x); @@ -2286,8 +2878,11 @@ int main(int argc, char *argv[]) int min_num_peaks_per_pixel = 3; int min_num_peaks_per_panel = 100; int only_best_distance = 0; + int enforce_cspad_layout = 0; int nostretch = 0; int individual_coffset = 0; + int no_cspad = 0; + int error_maps = 1; double max_peak_dist = 4.0; struct detector *det = NULL; @@ -2298,25 +2893,26 @@ int main(int argc, char *argv[]) const struct option longopts[] = { /* Options with long and short versions */ - {"help", 0, NULL, 'h'}, - {"version", 0, NULL, 10 }, - {"input", 1, NULL, 'i'}, - {"output", 1, NULL, 'o'}, - {"geometry", 1, NULL, 'g'}, - {"quadrants", 1, NULL, 'q'}, - {"connected", 1, NULL, 'c'}, - {"min-num-peaks-per-pixel",1, NULL, 'x'}, - {"min-num-peaks-per-panel",1, NULL, 'p'}, - {"most-few-clen", 0, NULL, 'l'}, - {"max-peak-dist", 1, NULL, 'm'}, - {"individual-dist-offset", 0, NULL, 's'}, - + {"help", 0, NULL, 'h'}, + {"version", 0, NULL, 10 }, + {"input", 1, NULL, 'i'}, + {"output", 1, NULL, 'o'}, + {"geometry", 1, NULL, 'g'}, + {"quadrants", 1, NULL, 'q'}, + {"connected", 1, NULL, 'c'}, + {"min-num-peaks-per-pixel", 1, NULL, 'x'}, + {"min-num-peaks-per-panel", 1, NULL, 'p'}, + {"most-few-clen", 0, NULL, 'l'}, + {"max-peak-dist", 1, NULL, 'm'}, + {"individual-dist-offset", 0, NULL, 's'}, /* Long-only options with no arguments */ - {"no-stretch", 0, &nostretch, 1}, - + {"no-stretch", 0, &nostretch, 1}, + {"no-error-maps", 0, &error_maps, 0}, + {"enforce-cspad-layout", 0, &enforce_cspad_layout, 1}, + {"no-cspad", 0, &no_cspad, 1}, - {0, 0, NULL, 0} + {0, 0, NULL, 0} }; /* Short options */ @@ -2379,6 +2975,7 @@ int main(int argc, char *argv[]) case 's' : individual_coffset = 1; break; + } } @@ -2398,7 +2995,8 @@ int main(int argc, char *argv[]) } if ( quadrant_coll_name == NULL ) { - ERROR("You must provide a rigid group collection for quadrants.\n"); + ERROR("You must provide a rigid group collection for " + "quadrants.\n"); return 1; } @@ -2417,10 +3015,11 @@ int main(int argc, char *argv[]) return 1; } - connected = find_rigid_group_collection_by_name(det, connected_coll_name); + connected = find_rigid_group_collection_by_name(det, + connected_coll_name); if ( connected == NULL ) { - ERROR("Cannot find rigid group collection for connected asics: %s\n", - connected_coll_name); + ERROR("Cannot find rigid group collection for connected " + "asics: %s\n", connected_coll_name); return 1; } @@ -2430,10 +3029,12 @@ int main(int argc, char *argv[]) strcat(command_line, buffer); } + g_type_init(); ret_val = optimize_geometry(infile, outfile, geometry_filename, det, quadrants, connected, min_num_peaks_per_pixel, min_num_peaks_per_panel, only_best_distance, - nostretch, individual_coffset, + nostretch, individual_coffset, error_maps, + enforce_cspad_layout, no_cspad, max_peak_dist, command_line); return ret_val; diff --git a/src/hdfsee-render.c b/src/hdfsee-render.c index abbb2ccd..71e2e1e8 100644 --- a/src/hdfsee-render.c +++ b/src/hdfsee-render.c @@ -46,16 +46,15 @@ #include <image.h> static float *get_binned_panel(struct image *image, int binning, - struct panel *p, double *max, int *pw, int *ph) + int pi, double *max, int *pw, int *ph) { float *data; int x, y; int w, h; int fw; - float *in; + struct panel *p = &image->det->panels[pi]; - fw = image->width; - in = image->data; + fw = p->max_fs - p->min_fs + 1; /* Some pixels might get discarded */ w = (p->max_fs - p->min_fs + 1) / binning; @@ -80,31 +79,13 @@ static float *get_binned_panel(struct image *image, int binning, double v; int fs, ss; - int tbad = 0; - fs = binning*x+xb+p->min_fs; - ss = binning*y+yb+p->min_ss; - v = in[fs+ss*fw]; + fs = binning*x+xb; + ss = binning*y+yb; + v = image->dp[pi][fs+ss*fw]; total += v; - if ( in_bad_region(image->det, fs, ss) ) tbad = 1; - - if ( image->flags != NULL ) { - - uint16_t flags = image->flags[fs+ss*fw]; - - if ( !((flags & image->det->mask_good) - == image->det->mask_good) ) { - tbad = 1; - } - - if ( flags & image->det->mask_bad ) { - tbad = 1; - } - - } - - if ( tbad ) bad = 1; + if ( image->bad[pi][fs+ss*fw] ) bad = 1; } } @@ -205,8 +186,7 @@ GdkPixbuf **render_panels(struct image *image, max = 0.0; for ( i=0; i<np; i++ ) { double this_max = 0.0; - hdrs[i] = get_binned_panel(image, binning, - &image->det->panels[i], &this_max, + hdrs[i] = get_binned_panel(image, binning, i, &this_max, &ws[i], &hs[i]); if ( this_max > max ) max = this_max; } @@ -246,7 +226,7 @@ GdkPixbuf *render_get_colour_scale(size_t w, size_t h, int scale) data = malloc(3*w*h); if ( data == NULL ) return NULL; - max = h; + max = h-(h/6); for ( y=0; y<h; y++ ) { diff --git a/src/hdfsee.c b/src/hdfsee.c index 89b8bf95..baf439bd 100644 --- a/src/hdfsee.c +++ b/src/hdfsee.c @@ -75,9 +75,11 @@ static void show_help(const char *s) " black-blue-pink-red-orange-\n" " -yellow-white.\n" " -e, --image=<element> Start up displaying this image from the\n" -" HDF5 file. Example: /data/data0.\n" -" (Only used when a geometry file is not" -" provided. See option -g)" +" HDF5 file. When this option is used,\n" +" information about the data layout\n" +" from the geometry file is ignored (See\n" +" manual page).\n" +" Example: /data/data0.\n" " --event=<event code> Event to show from multi-event file.\n" " -g, --geometry=<filename> Use geometry from file for display.\n" " (When this option is used, the value of\n" @@ -300,6 +302,11 @@ int main(int argc, char *argv[]) return 1; } + if ( event != NULL && geom_filename == NULL) { + ERROR("The '--event' option requires geometry file\n"); + return 1; + } + if ( cscale == NULL ) cscale = strdup("colour"); if ( strcmp(cscale, "mono") == 0 ) { colscale = SCALE_MONO; diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c deleted file mode 100644 index 7f15f54d..00000000 --- a/src/hrs-scaling.c +++ /dev/null @@ -1,559 +0,0 @@ -/* - * hrs-scaling.c - * - * Intensity scaling using generalised HRS target function - * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2010-2014 Thomas White <taw@physics.org> - * - * 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 <stdlib.h> -#include <assert.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_linalg.h> -#include <gsl/gsl_eigen.h> - -#include "image.h" -#include "peaks.h" -#include "symmetry.h" -#include "geometry.h" -#include "cell.h" -#include "utils.h" -#include "reflist.h" - - -/* Minimum partiality of a reflection for it to be used for scaling */ -#define MIN_PART_SCALE (0.05) - -/* Minimum partiality of a reflection for it to be merged */ -#define MIN_PART_MERGE (0.05) - -/* Maximum number of iterations of scaling per macrocycle. */ -#define MAX_CYCLES (10) - - -struct scale_queue_args -{ - RefList *reference; - Crystal **crystals; - int n_started; - PartialityModel pmodel; -}; - - -struct scale_worker_args -{ - Crystal *crystal; - RefList *reference; - PartialityModel pmodel; -}; - - -static void *create_scale_job(void *vqargs) -{ - struct scale_worker_args *wargs; - struct scale_queue_args *qargs = vqargs; - - wargs = malloc(sizeof(struct scale_worker_args)); - wargs->reference = qargs->reference; - wargs->pmodel = qargs->pmodel; - - wargs->crystal = qargs->crystals[qargs->n_started++]; - - return wargs; -} - - -static void run_scale_job(void *vwargs, int cookie) -{ - struct scale_worker_args *wargs = vwargs; - Crystal *cr = wargs->crystal; - RefList *reference = wargs->reference; - Reflection *refl; - RefListIterator *iter; - double num = 0.0; - double den = 0.0; - double g; - - for ( refl = first_refl(crystal_get_reflections(cr), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l; - double Ih, Ihl, corr; - Reflection *r; - - if ( (get_partiality(refl) < MIN_PART_SCALE) - || (get_intensity(refl) < 3.0*get_esd_intensity(refl)) ) { - continue; - } - - /* Look up by asymmetric indices */ - get_indices(refl, &h, &k, &l); - r = find_refl(reference, h, k, l); - if ( r == NULL ) continue; - - Ih = get_intensity(r); - - corr = get_partiality(refl) * get_lorentz(refl); - - Ihl = get_intensity(refl) / corr; - - num += Ih * Ihl; - den += Ih * Ih; - - } - - g = num / den; - crystal_set_osf(cr, g); /* If it's NaN, it'll get rejected later */ -} - - -static void finalise_scale_job(void *vqargs, void *vwargs) -{ - struct scale_worker_args *wargs = vwargs; - free(wargs); -} - - -static void iterate_scale(Crystal **crystals, int n, RefList *reference, - int n_threads, PartialityModel pmodel) -{ - struct scale_queue_args qargs; - - assert(reference != NULL); - - qargs.reference = reference; - qargs.n_started = 0; - qargs.crystals = crystals; - qargs.pmodel = pmodel; - - run_threads(n_threads, run_scale_job, create_scale_job, - finalise_scale_job, &qargs, n, 0, 0, 0); -} - - -struct merge_queue_args -{ - RefList *full; - pthread_rwlock_t full_lock; - Crystal **crystals; - int n_started; - PartialityModel pmodel; -}; - - -struct merge_worker_args -{ - Crystal *crystal; - RefList *full; - pthread_rwlock_t *full_lock; - PartialityModel pmodel; -}; - - -static void *create_merge_job(void *vqargs) -{ - struct merge_worker_args *wargs; - struct merge_queue_args *qargs = vqargs; - - wargs = malloc(sizeof(struct merge_worker_args)); - wargs->full = qargs->full; - wargs->full_lock = &qargs->full_lock; - wargs->pmodel = qargs->pmodel; - - wargs->crystal = qargs->crystals[qargs->n_started++]; - - return wargs; -} - - -static void run_merge_job(void *vwargs, int cookie) -{ - struct merge_worker_args *wargs = vwargs; - Crystal *cr = wargs->crystal; - RefList *full = wargs->full; - Reflection *refl; - RefListIterator *iter; - double G; - - /* If this crystal's scaling was dodgy, it doesn't contribute to the - * merged intensities */ - if ( crystal_get_user_flag(cr) != 0 ) return; - - G = crystal_get_osf(cr); - - for ( refl = first_refl(crystal_get_reflections(cr), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - Reflection *f; - signed int h, k, l; - double num, den; - int red; - double Ihl, corr; - - if ( get_partiality(refl) < MIN_PART_MERGE ) continue; - - get_indices(refl, &h, &k, &l); - pthread_rwlock_rdlock(wargs->full_lock); - f = find_refl(full, h, k, l); - if ( f == NULL ) { - - /* Swap read lock for write lock */ - pthread_rwlock_unlock(wargs->full_lock); - pthread_rwlock_wrlock(wargs->full_lock); - - /* In the gap between the unlock and the wrlock, the - * reflection might have been created by another thread. - * So, we must check again */ - f = find_refl(full, h, k, l); - if ( f == NULL ) { - f = add_refl(full, h, k, l); - lock_reflection(f); - pthread_rwlock_unlock(wargs->full_lock); - num = 0.0; - den = 0.0; - red = 0; - - } else { - - /* Someone else created it */ - lock_reflection(f); - pthread_rwlock_unlock(wargs->full_lock); - num = get_temp1(f); - den = get_temp2(f); - red = get_redundancy(f); - - } - - } else { - - lock_reflection(f); - pthread_rwlock_unlock(wargs->full_lock); - num = get_temp1(f); - den = get_temp2(f); - red = get_redundancy(f); - - } - - corr = get_partiality(refl) * get_lorentz(refl); - - Ihl = get_intensity(refl) / corr; - - num += Ihl / G; - den += 1.0; - red++; - - set_temp1(f, num); - set_temp2(f, den); - set_redundancy(f, red); - unlock_reflection(f); - } -} - - -static void finalise_merge_job(void *vqargs, void *vwargs) -{ - free(vwargs); -} - - -static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads, - PartialityModel pmodel) -{ - RefList *full; - struct merge_queue_args qargs; - Reflection *refl; - RefListIterator *iter; - - full = reflist_new(); - - qargs.full = full; - qargs.n_started = 0; - qargs.crystals = crystals; - qargs.pmodel = pmodel; - pthread_rwlock_init(&qargs.full_lock, NULL); - - run_threads(n_threads, run_merge_job, create_merge_job, - finalise_merge_job, &qargs, n, 0, 0, 0); - - pthread_rwlock_destroy(&qargs.full_lock); - - for ( refl = first_refl(full, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double Ih; - Ih = get_temp1(refl) / get_temp2(refl); - set_intensity(refl, Ih); - } - - return full; -} - - - -struct esd_queue_args -{ - RefList *full; - Crystal **crystals; - int n_started; - PartialityModel pmodel; -}; - - -struct esd_worker_args -{ - Crystal *crystal; - RefList *full; - PartialityModel pmodel; -}; - - -static void *create_esd_job(void *vqargs) -{ - struct esd_worker_args *wargs; - struct esd_queue_args *qargs = vqargs; - - wargs = malloc(sizeof(struct esd_worker_args)); - wargs->full = qargs->full; - wargs->pmodel = qargs->pmodel; - - wargs->crystal = qargs->crystals[qargs->n_started++]; - - return wargs; -} - - -static void run_esd_job(void *vwargs, int cookie) -{ - struct esd_worker_args *wargs = vwargs; - Crystal *cr = wargs->crystal; - RefList *full = wargs->full; - Reflection *refl; - RefListIterator *iter; - double G; - - /* If this crystal's scaling was dodgy, it doesn't contribute to the - * merged intensities */ - if ( crystal_get_user_flag(cr) != 0 ) return; - - G = crystal_get_osf(cr); - - for ( refl = first_refl(crystal_get_reflections(cr), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - Reflection *f; - signed int h, k, l; - double num; - double Ihl, Ih, corr; - - if ( get_partiality(refl) < MIN_PART_MERGE ) continue; - - get_indices(refl, &h, &k, &l); - f = find_refl(full, h, k, l); - assert(f != NULL); - - lock_reflection(f); - - num = get_temp1(f); - - corr = get_partiality(refl) * get_lorentz(refl); - - Ih = get_intensity(f); - Ihl = get_intensity(refl) / (G*corr); - - num += pow(Ihl - Ih, 2.0); - - set_temp1(f, num); - unlock_reflection(f); - } -} - - -static void finalise_esd_job(void *vqargs, void *vwargs) -{ - free(vwargs); -} - - -static void calculate_esds(Crystal **crystals, int n, RefList *full, - int n_threads, int min_red, PartialityModel pmodel) -{ - struct esd_queue_args qargs; - Reflection *refl; - RefListIterator *iter; - - qargs.full = full; - qargs.n_started = 0; - qargs.crystals = crystals; - qargs.pmodel = pmodel; - - for ( refl = first_refl(full, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - set_temp1(refl, 0.0); - set_temp2(refl, 0.0); - } - - run_threads(n_threads, run_esd_job, create_esd_job, - finalise_esd_job, &qargs, n, 0, 0, 0); - - for ( refl = first_refl(full, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double esd; - int red = get_redundancy(refl); - esd = sqrt(get_temp1(refl)); - esd /= (double)red; - set_esd_intensity(refl, esd); - - if ( red < min_red ) { - set_redundancy(refl, 0); - } - } -} - - -static void reject_outliers(double *old_osfs, int n, Crystal **crystals) -{ - int i; - - for ( i=0; i<n; i++ ) { - double osf = crystal_get_osf(crystals[i]); - if ( isnan(osf) || (osf < 0.0) || (osf > 3.0) ) { - crystal_set_user_flag(crystals[i], 1); - } - } -} - - -static int test_convergence(double *old_osfs, int n, Crystal **crystals) -{ - int i; - double total_change = 0.0; - double mean_change; - int n_change = 0; - - for ( i=0; i<n; i++ ) { - if ( crystal_get_user_flag(crystals[i]) == 0 ) { - double new_osf = crystal_get_osf(crystals[i]); - total_change += fabs(new_osf - old_osfs[i]); - n_change++; - } - } - mean_change = total_change / n_change; - - STATUS("Mean OSF change = %f\n", mean_change); - - return mean_change < 0.01; -} - - -/* Scale the stack of images */ -RefList *scale_intensities(Crystal **crystals, int n, - int n_threads, int noscale, PartialityModel pmodel, - int min_redundancy) -{ - int i; - RefList *full = NULL; - double *old_osfs; - int done; - - for ( i=0; i<n; i++ ) { - crystal_set_user_flag(crystals[i], 0); - crystal_set_osf(crystals[i], 1.0); - } - - if ( noscale ) { - full = lsq_intensities(crystals, n, n_threads, pmodel); - calculate_esds(crystals, n, full, n_threads, min_redundancy, - pmodel); - return full; - } - - /* Create an initial list to refine against */ - full = lsq_intensities(crystals, n, n_threads, pmodel); - - old_osfs = malloc(n*sizeof(double)); - if ( old_osfs == NULL ) return NULL; - - /* Iterate */ - i = 0; - do { - - double total_sf = 0.0; - int n_sf = 0; - double norm_sf; - int j; - - for ( j=0; j<n; j++ ) { - old_osfs[j] = crystal_get_osf(crystals[j]); - crystal_set_user_flag(crystals[j], 0); - } - - iterate_scale(crystals, n, full, n_threads, pmodel); - - /* Normalise the scale factors */ - for ( j=0; j<n; j++ ) { - double osf = crystal_get_osf(crystals[j]); - if ( !isnan(osf) ) { - total_sf += osf; - n_sf++; - } - } - norm_sf = total_sf / n_sf; - for ( j=0; j<n; j++ ) { - crystal_set_osf(crystals[j], - crystal_get_osf(crystals[j])/norm_sf); - } - - reject_outliers(old_osfs, n, crystals); - done = test_convergence(old_osfs, n, crystals); - - /* Generate list for next iteration */ - reflist_free(full); - full = lsq_intensities(crystals, n, n_threads, pmodel); - - i++; - - } while ( !done && (i < MAX_CYCLES) ); - - if ( i == MAX_CYCLES ) { - ERROR("WARNING: Scaling did not converge.\n"); - } - - calculate_esds(crystals, n, full, n_threads, min_redundancy, pmodel); - - free(old_osfs); - return full; -} diff --git a/src/im-sandbox.c b/src/im-sandbox.c index ba709f6f..df50199c 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -3,13 +3,13 @@ * * Sandbox for indexing * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2014 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * 2014 Valerio Mariani * 2011 Richard Kirian * 2012 Lorenzo Galli @@ -48,6 +48,7 @@ #include <signal.h> #include <sys/stat.h> #include <assert.h> +#include <sys/mman.h> #ifdef HAVE_CLOCK_GETTIME #include <time.h> @@ -85,6 +86,12 @@ struct sb_reader }; +struct sb_shm +{ + pthread_mutex_t term_lock; +}; + + struct sandbox { pthread_mutex_t lock; @@ -96,7 +103,6 @@ struct sandbox int n_hadcrystals_last_stats; int n_crystals_last_stats; int t_last_stats; - int suspend_stats; struct index_args *iargs; @@ -110,6 +116,8 @@ struct sandbox struct filename_plus_event **last_filename; int serial; + struct sb_shm *shared; + char *tmpdir; struct sb_reader *reader; @@ -359,7 +367,7 @@ static int read_fpe_data(struct buffer_data *bd) static void run_work(const struct index_args *iargs, int filename_pipe, int results_pipe, Stream *st, - int cookie, const char *tmpdir) + int cookie, const char *tmpdir, pthread_mutex_t *term_lock) { FILE *fh; int allDone = 0; @@ -494,7 +502,7 @@ static void run_work(const struct index_args *iargs, pargs.n_crystals = 0; process_image(iargs, &pargs, st, cookie, tmpdir, - results_pipe, ser); + results_pipe, ser, term_lock); /* Request another image */ c = sprintf(buf, "%i\n", pargs.n_crystals); @@ -815,11 +823,12 @@ static void start_worker_process(struct sandbox *sb, int slot) st = open_stream_fd_for_write(stream_pipe[1]); run_work(sb->iargs, filename_pipe[0], result_pipe[1], - st, slot, tmp); + st, slot, tmp, &sb->shared->term_lock); close_stream(st); //close(filename_pipe[0]); close(result_pipe[1]); + munmap(sb->shared, sizeof(struct sb_shm)); free(sb); @@ -894,6 +903,39 @@ static void handle_zombie(struct sandbox *sb) } +static int setup_shm(struct sandbox *sb) +{ + pthread_mutexattr_t attr; + + sb->shared = mmap(NULL, sizeof(struct sb_shm), PROT_READ | PROT_WRITE, + MAP_SHARED | MAP_ANONYMOUS, -1, 0); + + if ( sb->shared == MAP_FAILED ) { + ERROR("SHM setup failed: %s\n", strerror(errno)); + return 1; + } + + if ( pthread_mutexattr_init(&attr) ) { + ERROR("Failed to initialise mutex attr.\n"); + return 1; + } + + if ( pthread_mutexattr_setpshared(&attr, PTHREAD_PROCESS_SHARED) ) { + ERROR("Failed to set process shared attribute.\n"); + return 1; + } + + if ( pthread_mutex_init(&sb->shared->term_lock, &attr) ) { + ERROR("Terminal lock setup failed.\n"); + return 1; + } + + pthread_mutexattr_destroy(&attr); + + return 0; +} + + void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, Stream *stream, const char *tempdir) @@ -930,7 +972,6 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->n_hadcrystals_last_stats = 0; sb->n_crystals_last_stats = 0; sb->t_last_stats = get_monotonic_seconds(); - sb->suspend_stats = 0; sb->n_proc = n_proc; sb->iargs = iargs; sb->serial = 1; @@ -939,6 +980,12 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->reader->fhs = NULL; sb->reader->stream = stream; + if ( setup_shm(sb) ) { + ERROR("Failed to set up SHM.\n"); + free(sb); + return; + } + sb->stream_pipe_write = calloc(n_proc, sizeof(int)); if ( sb->stream_pipe_write == NULL ) { ERROR("Couldn't allocate memory for pipes.\n"); @@ -1106,34 +1153,20 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, chomp(results); - if ( strcmp(results, "SUSPEND") == 0 ) { - sb->suspend_stats++; - continue; /* Do not send next filename */ - } else if ( strcmp(results, "RELEASE") == 0 ) { - if ( sb->suspend_stats > 0 ) { - sb->suspend_stats--; - } else { - ERROR("RELEASE before SUSPEND.\n"); + strtol(results, &eptr, 10); + if ( eptr == results ) { + if ( strlen(results) > 0 ) { + ERROR("Invalid result '%s'\n", + results); } - continue; /* Do not send next filename */ } else { - strtol(results, &eptr, 10); - if ( eptr == results ) { - if ( strlen(results) > 0 ) { - ERROR("Invalid result '%s'\n", - results); - } - } else { - - int nc = atoi(results); - sb->n_crystals += nc; - if ( nc > 0 ) { - sb->n_hadcrystals++; - } - sb->n_processed++; - + int nc = atoi(results); + sb->n_crystals += nc; + if ( nc > 0 ) { + sb->n_hadcrystals++; } + sb->n_processed++; } @@ -1210,8 +1243,8 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, /* Update progress */ lock_sandbox(sb); tNow = get_monotonic_seconds(); - if ( !sb->suspend_stats - && (tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS) ) + r = pthread_mutex_trylock(&sb->shared->term_lock); + if ((r==0) && (tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS)) { STATUS("%4i indexable out of %4i processed (%4.1f%%), " @@ -1228,7 +1261,9 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->n_crystals_last_stats = sb->n_crystals; sb->t_last_stats = tNow; + } + if ( r == 0 ) pthread_mutex_unlock(&sb->shared->term_lock); unlock_sandbox(sb); allDone = 1; @@ -1264,6 +1299,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, free(sb->result_fhs); free(sb->pids); free(sb->tmpdir); + munmap(sb->shared, sizeof(struct sb_shm)); pthread_mutex_destroy(&sb->lock); diff --git a/src/indexamajig.c b/src/indexamajig.c index 44f067e2..c7e4a270 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -3,13 +3,13 @@ * * Index patterns, output hkl+intensity etc. * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2014 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * 2011 Richard Kirian * 2012 Lorenzo Galli * 2012 Chunhong Yoon @@ -97,11 +97,9 @@ static void show_help(const char *s) " zaef : Use Zaefferer (2000) gradient detection.\n" " This is the default method.\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" " Default: /processing/hitfinder/peakinfo\n" -" --cxi-hdf5-peaks Peaks in the HDF5 file are in CXI file format.\n" -" Only used in conjunction with the --hdf5-peaks,\n" -" ignored otherwise." " --integration=<meth> Perform final pattern integration using <meth>.\n" "\n\n" "For more control over the process, you might need:\n\n" @@ -127,6 +125,10 @@ static void show_help(const char *s) " --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" @@ -145,6 +147,7 @@ static void show_help(const char *s) " --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" ); } @@ -214,8 +217,7 @@ int main(int argc, char *argv[]) iargs.det = NULL; iargs.peaks = PEAK_ZAEF; iargs.beam = &beam; - iargs.hdf5_peak_path = strdup("/processing/hitfinder/peakinfo"); - iargs.cxi_hdf5_peaks = 0; + iargs.hdf5_peak_path = NULL; iargs.copyme = NULL; iargs.pk_inn = -1.0; iargs.pk_mid = -1.0; @@ -238,6 +240,10 @@ int main(int argc, char *argv[]) iargs.int_meth = integration_method("rings-nocen", NULL); iargs.push_res = 0.0; iargs.highres = +INFINITY; + iargs.fix_profile_r = -1.0; + iargs.fix_bandwidth = -1.0; + iargs.fix_divergence = -1.0; + iargs.predict_refine = 1; /* Long options */ const struct option longopts[] = { @@ -264,7 +270,7 @@ int main(int argc, char *argv[]) {"no-use-saturated", 0, &iargs.use_saturated, 0}, {"no-revalidate", 0, &iargs.no_revalidate, 1}, {"check-hdf5-snr", 0, &iargs.check_hdf5_snr, 1}, - {"cxi-hdf5-peaks", 0, &iargs.cxi_hdf5_peaks, 1}, + {"no-refine", 0, &iargs.predict_refine, 0}, /* Long-only options which don't actually do anything */ {"no-sat-corr", 0, &iargs.satcorr, 0}, @@ -293,6 +299,9 @@ int main(int argc, char *argv[]) {"res-push", 1, NULL, 19}, /* compat */ {"peak-radius", 1, NULL, 20}, {"highres", 1, NULL, 21}, + {"fix-profile-radius", 1, NULL, 22}, + {"fix-bandwidth", 1, NULL, 23}, + {"fix-divergence", 1, NULL, 24}, {0, 0, NULL, 0} }; @@ -440,6 +449,28 @@ int main(int argc, char *argv[]) iargs.highres = 1.0 / (iargs.highres/1e10); break; + case 22 : + if ( sscanf(optarg, "%f", &iargs.fix_profile_r) != 1 ) { + ERROR("Invalid value for " + "--fix-profile-radius\n"); + return 1; + } + break; + + case 23 : + if ( sscanf(optarg, "%f", &iargs.fix_bandwidth) != 1 ) { + ERROR("Invalid value for --fix-bandwidth\n"); + return 1; + } + break; + + case 24 : + if ( sscanf(optarg, "%f", &iargs.fix_divergence) != 1 ) { + ERROR("Invalid value for --fix-divergence\n"); + return 1; + } + break; + case 0 : break; @@ -481,12 +512,23 @@ int main(int argc, char *argv[]) iargs.peaks = PEAK_ZAEF; } else if ( strcmp(speaks, "hdf5") == 0 ) { iargs.peaks = PEAK_HDF5; + } else if ( strcmp(speaks, "cxi") == 0 ) { + iargs.peaks = PEAK_CXI; } else { ERROR("Unrecognised peak detection method '%s'\n", speaks); return 1; } free(speaks); + /* Set default path for peaks, if appropriate */ + if ( iargs.hdf5_peak_path == NULL ) { + if ( iargs.peaks == PEAK_HDF5 ) { + iargs.hdf5_peak_path = strdup("/processing/hitfinder/peakinfo"); + } else if ( iargs.peaks == PEAK_CXI ) { + iargs.hdf5_peak_path = strdup("/entry_1/result_1"); + } + } + if ( prefix == NULL ) { prefix = strdup(""); } else { @@ -653,11 +695,6 @@ int main(int argc, char *argv[]) free(int_diag); - if ( (n_proc > 1) && (iargs.int_diag != INTDIAG_NONE) ) { - n_proc = 1; - STATUS("Ignored \"-j\" because you used --int-diag.\n"); - } - } st = open_stream_for_write_2(outfile, geom_filename, argc, argv); diff --git a/src/list_events.c b/src/list_events.c new file mode 100644 index 00000000..1ef82701 --- /dev/null +++ b/src/list_events.c @@ -0,0 +1,195 @@ +/* + * list_events.c + * + * Generate event lists + * + * Copyright © 2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2015 Thomas White <taw@physics.org> + * + * 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 <stdarg.h> +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <unistd.h> +#include <getopt.h> + +#include "version.h" +#include "utils.h" +#include "detector.h" +#include "hdf5-file.h" + + +static void show_help(const char *s) +{ + printf("Syntax: %s [options] -i files.lst -o events.lst " + "-g geometry.geom\n\n", s); + printf( +"Generate event lists.\n" +"\n" +" -h, --help Display this help message.\n" +" --version Print CrystFEL version number and exit.\n" +"\n" +" -i, --input=<file> Input filename (list of multi-event filenames).\n" +" -g, --geometry=<file> Get data layout from geometry file.\n" +" -o, --output=<file> Output filename (list of events).\n" +); +} + + +int main(int argc, char *argv[]) +{ + int c; + char *input = NULL; + char *output = NULL; + char *geom = NULL; + char *rval; + FILE *ifh; + FILE *ofh; + struct detector *det; + + /* Long options */ + const struct option longopts[] = { + {"help", 0, NULL, 'h'}, + {"version", 0, NULL, 2 }, + {"input", 1, NULL, 'i'}, + {"geometry", 1, NULL, 'g'}, + {"output", 1, NULL, 'o'}, + {0, 0, NULL, 0} + }; + + /* Short options */ + while ((c = getopt_long(argc, argv, "hi:g:o:", + longopts, NULL)) != -1) { + + switch (c) { + + case 'h' : + show_help(argv[0]); + return 0; + + case 2 : + printf("CrystFEL: " CRYSTFEL_VERSIONSTRING "\n"); + printf(CRYSTFEL_BOILERPLATE"\n"); + return 0; + + case 'o' : + output = strdup(optarg); + break; + + case 'i' : + input = strdup(optarg); + break; + + case 'g' : + geom = strdup(optarg); + break; + + case 0 : + break; + + case '?' : + break; + + default : + ERROR("Unhandled option '%c'\n", c); + break; + + } + + } + + if ( (input == NULL) || (output == NULL) || (geom == NULL) ) { + ERROR("You must specify at least the input, output and geometry" + " filenames.\n"); + return 1; + } + + ifh = fopen(input, "r"); + if ( ifh == NULL ) { + ERROR("Couldn't open '%s'\n", input); + return 1; + } + + ofh = fopen(output, "w"); + if ( ofh == NULL ) { + ERROR("Couldn't open '%s'\n", output); + return 1; + } + + det = get_detector_geometry(geom, NULL); + if ( det == NULL ) { + ERROR("Failed to read '%s'\n", geom); + return 1; + } + + if ( (det->path_dim == 0) && (det->dim_dim == 0) ) { + ERROR("This does not look like a multi-event geometry file.\n"); + ERROR("Are you sure you need to use list_events instead of " + "just 'find' or 'ls'?\n"); + return 1; + } + + do { + + char filename[1024]; + int i; + + rval = fgets(filename, 1024, ifh); + if ( rval != NULL ) { + + struct event_list *evlist; + struct hdfile *hdfile; + + chomp(filename); + + hdfile = hdfile_open(filename); + if ( hdfile == NULL ) { + ERROR("Failed to open '%s'\n", filename); + ERROR("Aborting creation of event list.\n"); + return 1; + } + + evlist = fill_event_list(hdfile, det); + + for ( i=0; i<evlist->num_events; i++ ) { + char *str = get_event_string(evlist->events[i]); + fprintf(ofh, "%s %s\n", filename, str); + free(str); + } + + STATUS("%i events found in %s\n", evlist->num_events, + filename); + + free_event_list(evlist); + hdfile_close(hdfile); + } + + } while ( rval != NULL ); + + return 0; +} diff --git a/src/merge.c b/src/merge.c new file mode 100644 index 00000000..1595aea3 --- /dev/null +++ b/src/merge.c @@ -0,0 +1,248 @@ +/* + * merge.c + * + * Parallel weighted merging of intensities + * + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2015 Thomas White <taw@physics.org> + * + * 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 <stdlib.h> +#include <assert.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_eigen.h> +#include <gsl/gsl_fit.h> + +#include "image.h" +#include "peaks.h" +#include "symmetry.h" +#include "geometry.h" +#include "cell.h" +#include "utils.h" +#include "reflist.h" +#include "cell-utils.h" + + +/* Minimum partiality of a reflection for it to be merged */ +#define MIN_PART_MERGE (0.05) + + +struct merge_queue_args +{ + RefList *full; + pthread_rwlock_t full_lock; + Crystal **crystals; + int n_started; + PartialityModel pmodel; + double push_res; +}; + + +struct merge_worker_args +{ + struct merge_queue_args *qargs; + Crystal *crystal; + int crystal_number; +}; + + +static void *create_merge_job(void *vqargs) +{ + struct merge_worker_args *wargs; + struct merge_queue_args *qargs = vqargs; + + wargs = malloc(sizeof(struct merge_worker_args)); + wargs->qargs = qargs; + wargs->crystal_number = qargs->n_started; + wargs->crystal = qargs->crystals[qargs->n_started++]; + + return wargs; +} + + +static void run_merge_job(void *vwargs, int cookie) +{ + struct merge_worker_args *wargs = vwargs; + Crystal *cr = wargs->crystal; + RefList *full = wargs->qargs->full; + double push_res = wargs->qargs->push_res; + Reflection *refl; + RefListIterator *iter; + double G, B; + + /* If this crystal's scaling was dodgy, it doesn't contribute to the + * merged intensities */ + if ( crystal_get_user_flag(cr) != 0 ) return; + + G = crystal_get_osf(cr); + B = crystal_get_Bfac(cr); + + for ( refl = first_refl(crystal_get_reflections(cr), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + Reflection *f; + signed int h, k, l; + double mean, sumweight, M2, temp, delta, R; + double corr, res, w, esd; + + if ( get_partiality(refl) < MIN_PART_MERGE ) continue; + + get_indices(refl, &h, &k, &l); + pthread_rwlock_rdlock(&wargs->qargs->full_lock); + f = find_refl(full, h, k, l); + if ( f == NULL ) { + + /* Swap read lock for write lock */ + pthread_rwlock_unlock(&wargs->qargs->full_lock); + pthread_rwlock_wrlock(&wargs->qargs->full_lock); + + /* In the gap between the unlock and the wrlock, the + * reflection might have been created by another thread. + * So, we must check again */ + f = find_refl(full, h, k, l); + if ( f == NULL ) { + f = add_refl(full, h, k, l); + lock_reflection(f); + pthread_rwlock_unlock(&wargs->qargs->full_lock); + set_intensity(f, 0.0); + set_temp1(f, 0.0); + set_temp2(f, 0.0); + + } else { + + /* Someone else created it */ + lock_reflection(f); + pthread_rwlock_unlock(&wargs->qargs->full_lock); + + } + + } else { + + lock_reflection(f); + pthread_rwlock_unlock(&wargs->qargs->full_lock); + + } + + mean = get_intensity(f); + sumweight = get_temp1(f); + M2 = get_temp2(f); + + res = resolution(crystal_get_cell(cr), h, k, l); + + if ( 2.0*res > crystal_get_resolution_limit(cr)+push_res ) { + unlock_reflection(f); + continue; + } + + /* Total (multiplicative) correction factor */ + corr = exp(-G) * exp(B*res*res) * get_lorentz(refl) + / get_partiality(refl); + if ( isnan(corr) ) { + ERROR("NaN corr:\n"); + ERROR("G = %f, B = %e\n", G, B); + ERROR("res = %e\n", res); + ERROR("p = %f\n", get_partiality(refl)); + } + + esd = get_esd_intensity(refl) * corr; + w = 1.0; + + /* Running mean and variance calculation */ + temp = w + sumweight; + delta = get_intensity(refl)*corr - mean; + R = delta * w / temp; + set_intensity(f, mean + R); + set_temp2(f, M2 + sumweight * delta * R); + set_temp1(f, temp); + set_redundancy(f, get_redundancy(f)+1); + unlock_reflection(f); + } +} + + +static void finalise_merge_job(void *vqargs, void *vwargs) +{ + free(vwargs); +} + + +RefList *merge_intensities(Crystal **crystals, int n, int n_threads, + PartialityModel pmodel, int min_meas, + double push_res) +{ + RefList *full; + RefList *full2; + struct merge_queue_args qargs; + Reflection *refl; + RefListIterator *iter; + + full = reflist_new(); + + qargs.full = full; + qargs.n_started = 0; + qargs.crystals = crystals; + qargs.pmodel = pmodel; + qargs.push_res = push_res; + pthread_rwlock_init(&qargs.full_lock, NULL); + + run_threads(n_threads, run_merge_job, create_merge_job, + finalise_merge_job, &qargs, n, 0, 0, 0); + + pthread_rwlock_destroy(&qargs.full_lock); + + /* Calculate ESDs from variances, including only reflections with + * enough measurements */ + full2 = reflist_new(); + if ( full2 == NULL ) return NULL; + for ( refl = first_refl(full, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double var; + int red; + + red = get_redundancy(refl); + var = get_temp2(refl) / get_temp1(refl); + set_esd_intensity(refl, sqrt(var)/sqrt(red)); + + if ( red >= min_meas ) { + + signed int h, k, l; + Reflection *r2; + + get_indices(refl, &h, &k, &l); + r2 = add_refl(full2, h, k, l); + copy_data(r2, refl); + } + } + + reflist_free(full); + return full2; +} diff --git a/src/hrs-scaling.h b/src/merge.h index 16368b79..b8af63b8 100644 --- a/src/hrs-scaling.h +++ b/src/merge.h @@ -1,13 +1,13 @@ /* - * hrs-scaling.h + * merge.h * - * Intensity scaling using generalised HRS target function + * Parallel weighted merging of intensities * - * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2013 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -26,8 +26,8 @@ * */ -#ifndef HRS_SCALING_H -#define HRS_SCALING_H +#ifndef MERGE_H +#define MERGE_H #ifdef HAVE_CONFIG_H @@ -39,9 +39,8 @@ #include "reflist.h" #include "geometry.h" -extern RefList *scale_intensities(Crystal **crystals, int n, int n_threads, - int noscale, PartialityModel pmodel, - int min_redundancy); +extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads, + PartialityModel pmodel, int min_meas, + double push_res); - -#endif /* HRS_SCALING_H */ +#endif /* MERGE */ diff --git a/src/partial_sim.c b/src/partial_sim.c index c0d9aaaa..867183da 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -241,6 +241,7 @@ static void show_help(const char *s) " --profile-radius Reciprocal space reflection profile radius in m^-1.\n" " Default 0.001e9 m^-1\n" " --photon-energy Photon energy in eV. Default 9000.\n" +" --really-random Be non-deterministic.\n" "\n" ); } @@ -414,6 +415,27 @@ static void finalise_job(void *vqargs, void *vwargs) } +static void fixup_geom(struct detector *det) +{ + int i; + + for ( i=0; i<det->n_panels; i++ ) { + det->panels[i].clen += det->panels[i].coffset; + } +} + + +static int geom_contains_references(struct detector *det) +{ + int i; + + for ( i=0; i<det->n_panels; i++ ) { + if ( det->panels[i].clen_from != NULL ) return 1; + } + return 0; +} + + int main(int argc, char *argv[]) { int c; @@ -707,6 +729,13 @@ int main(int argc, char *argv[]) ERROR("The value given on the command line " "(with --photon-energy) will be used instead.\n"); } + if ( geom_contains_references(det) ) { + ERROR("Geometry file contains a reference to an HDF5 location" + " for the camera length. Change it to a numerical value " + " and try again.\n"); + return 1; + } + fixup_geom(det); if ( sym_str == NULL ) sym_str = strdup("1"); sym = get_pointgroup(sym_str); @@ -762,6 +791,7 @@ int main(int argc, char *argv[]) free(output_file); image.det = det; + image.beam = &beam; image.width = det->max_fs + 1; image.height = det->max_ss + 1; @@ -776,6 +806,8 @@ int main(int argc, char *argv[]) image.num_peaks = 0; image.num_saturated_peaks = 0; image.spectrum_size = 0; + image.spectrum = NULL; + image.serial = 0; image.event = NULL; STATUS("Simulation parameters:\n"); diff --git a/src/partialator.c b/src/partialator.c index e5f965ed..9156536c 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -3,11 +3,11 @@ * * Scaling and post refinement for coherent nanocrystallography * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2013 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -50,11 +50,13 @@ #include <thread-pool.h> #include <reflist.h> #include <reflist-utils.h> +#include <cell.h> +#include <cell-utils.h> #include "version.h" #include "post-refinement.h" -#include "hrs-scaling.h" -#include "scaling-report.h" +#include "merge.h" +#include "rejection.h" static void show_help(const char *s) @@ -75,6 +77,7 @@ static void show_help(const char *s) " --min-measurements=<n> Minimum number of measurements to require.\n" " --no-polarisation Disable polarisation correction.\n" " --max-adu=<n> Saturation value of detector.\n" +" --push-res=<n> Merge higher than apparent resolution cutoff.\n" " -j <n> Run <n> analyses in parallel.\n"); } @@ -84,6 +87,7 @@ struct refine_args RefList *full; Crystal *crystal; PartialityModel pmodel; + int no_scale; struct prdata prdata; }; @@ -94,8 +98,11 @@ struct queue_args int n_done; Crystal **crystals; int n_crystals; - struct srdata *srdata; struct refine_args task_defaults; + double initial_residual; + double initial_free_residual; + double final_residual; + double final_free_residual; }; @@ -104,7 +111,8 @@ static void refine_image(void *task, int id) struct refine_args *pargs = task; Crystal *cr = pargs->crystal; - pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel); + pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel, + pargs->no_scale); } @@ -130,10 +138,10 @@ static void done_image(void *vqargs, void *task) struct refine_args *pargs = task; qargs->n_done++; - if ( pargs->prdata.refined ) { - qargs->srdata->n_refined += pargs->prdata.refined; - qargs->srdata->n_filtered += pargs->prdata.n_filtered; - } + qargs->initial_residual += pargs->prdata.initial_residual; + qargs->initial_free_residual += pargs->prdata.initial_free_residual; + qargs->final_residual += pargs->prdata.final_residual; + qargs->final_free_residual += pargs->prdata.final_free_residual; progress_bar(qargs->n_done, qargs->n_crystals, "Refining"); free(task); @@ -142,27 +150,29 @@ static void done_image(void *vqargs, void *task) static void refine_all(Crystal **crystals, int n_crystals, RefList *full, int nthreads, PartialityModel pmodel, - struct srdata *srdata) + int no_scale, + double *initial_residual, double *initial_free_residual, + double *final_residual, double *final_free_residual) { struct refine_args task_defaults; struct queue_args qargs; - /* If the partiality model is "p=1", this refinement is really, really - * easy... */ - if ( pmodel == PMODEL_UNITY ) return; - task_defaults.full = full; task_defaults.crystal = NULL; task_defaults.pmodel = pmodel; task_defaults.prdata.refined = 0; task_defaults.prdata.n_filtered = 0; + task_defaults.no_scale = no_scale; qargs.task_defaults = task_defaults; qargs.n_started = 0; qargs.n_done = 0; qargs.n_crystals = n_crystals; qargs.crystals = crystals; - qargs.srdata = srdata; + qargs.initial_residual = 0.0; + qargs.initial_free_residual = 0.0; + qargs.final_residual = 0.0; + qargs.final_free_residual = 0.0; /* Don't have threads which are doing nothing */ if ( n_crystals < nthreads ) nthreads = n_crystals; @@ -170,9 +180,10 @@ static void refine_all(Crystal **crystals, int n_crystals, run_threads(nthreads, refine_image, get_image, done_image, &qargs, n_crystals, 0, 0, 0); - STATUS("%5.2f eigenvalues filtered on final iteration per successfully " - "refined crystal\n", - (double)srdata->n_filtered/srdata->n_refined); + *initial_residual = qargs.initial_residual; + *initial_free_residual = qargs.initial_free_residual; + *final_residual = qargs.final_residual; + *final_free_residual = qargs.final_free_residual; } @@ -192,11 +203,29 @@ static void display_progress(int n_images, int n_crystals) static const char *str_flags(Crystal *cr) { - if ( crystal_get_user_flag(cr) ) { - return "N"; - } + switch ( crystal_get_user_flag(cr) ) { + + case 0 : + return "OK"; + + case 1 : + return "bad scaling"; + + case 2 : + return "not enough reflections"; + + case 3 : + return "PR solve failed"; + + case 4 : + return "PR lost too many reflections"; + + case 5 : + return "Early rejection"; - return "-"; + default : + return "Unknown flag"; + } } @@ -225,6 +254,186 @@ static RefList *apply_max_adu(RefList *list, double max_adu) } +static void skip_to_end(FILE *fh) +{ + int c; + do { + c = fgetc(fh); + } while ( (c != '\n') && (c != EOF) ); +} + + +static int set_initial_params(Crystal *cr, FILE *fh) +{ + if ( fh != NULL ) { + + int err; + int n; + float osf, B; + + err = fscanf(fh, "%i %f %f", &n, &osf, &B); + if ( err != 3 ) { + ERROR("Failed to read parameters.\n"); + return 1; + } + + crystal_set_osf(cr, osf); + crystal_set_Bfac(cr, B*1e-20); + + skip_to_end(fh); + + } else { + + crystal_set_osf(cr, 0.0); + crystal_set_Bfac(cr, 0.0); + + } + + return 0; +} + + +static void show_duds(Crystal **crystals, int n_crystals) +{ + int j; + int n_dud = 0; + int n_noscale = 0; + int n_noref = 0; + int n_solve = 0; + int n_early = 0; + int n_cc = 0; + + for ( j=0; j<n_crystals; j++ ) { + int flag; + flag = crystal_get_user_flag(crystals[j]); + if ( flag != 0 ) n_dud++; + switch ( flag ) { + + case 0: + break; + + case 1: + n_noscale++; + break; + + case 2: + n_noref++; + break; + + case 3: + n_solve++; + break; + + case 5: + n_early++; + break; + + case 6: + n_cc++; + break; + + default: + STATUS("Unknown flag %i\n", flag); + break; + } + } + + if ( n_dud ) { + STATUS("%i bad crystals:\n", n_dud); + STATUS(" %i scaling failed.\n", n_noscale); + STATUS(" %i not enough reflections.\n", n_noref); + STATUS(" %i solve failed.\n", n_solve); + STATUS(" %i early rejection.\n", n_early); + STATUS(" %i bad CC.\n", n_cc); + } +} + + +/* Flag a random 5% of reflections */ +static void select_free_reflections(RefList *list, gsl_rng *rng) +{ + Reflection *refl; + RefListIterator *iter; + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + set_flag(refl, random_flat(rng, 1.0) > 0.95); + } +} + + +static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr, + int fr) +{ + Reflection *refl; + RefListIterator *iter; + double G = crystal_get_osf(cr); + double B = crystal_get_Bfac(cr); + UnitCell *cell = crystal_get_cell(cr); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + double pobs, pcalc; + double res, corr, Ipart; + Reflection *match; + + if ( !get_flag(refl) ) continue; /* Not free-flagged */ + + get_indices(refl, &h, &k, &l); + res = resolution(cell, h, k, l); + if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue; + + match = find_refl(full, h, k, l); + if ( match == NULL ) continue; + + /* Calculated partiality */ + pcalc = get_partiality(refl); + + /* Observed partiality */ + corr = exp(-G) * exp(B*res*res) * get_lorentz(refl); + Ipart = get_intensity(refl) * corr; + pobs = Ipart / get_intensity(match); + + fprintf(fh, "%5i %4i %4i %4i %8.4f %8.3f %8.3f\n", + fr, h, k, l, 2*res/1e9, pcalc, pobs); + + } +} + + +static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals, + int iter) +{ + FILE *fh; + char tmp[256]; + int i; + + snprintf(tmp, 256, "pgraph-iter%i.dat", iter); + + fh = fopen(tmp, "w"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", tmp); + return; + } + + fprintf(fh, " fr h k l 1/d(nm) pcalc pobs\n"); + + for ( i=0; i<n_crystals; i++ ) { + if ( crystal_get_user_flag(crystals[i]) != 0 ) continue; + write_to_pgraph(fh, crystal_get_reflections(crystals[i]), full, + crystals[i], i); + } + + fclose(fh); + +} + + int main(int argc, char *argv[]) { int c; @@ -240,17 +449,19 @@ int main(int argc, char *argv[]) int n_images = 0; int n_crystals = 0; char cmdline[1024]; - SRContext *sr; - int noscale = 0; + int no_scale = 0; Stream *st; Crystal **crystals; char *pmodel_str = NULL; PartialityModel pmodel = PMODEL_SCSPHERE; int min_measurements = 2; char *rval; - struct srdata srdata; int polarisation = 1; double max_adu = +INFINITY; + char *sparams_fn = NULL; + FILE *sparams_fh; + double push_res = 0.0; + gsl_rng *rng; /* Long options */ const struct option longopts[] = { @@ -266,8 +477,11 @@ int main(int argc, char *argv[]) {"min-measurements", 1, NULL, 2}, {"max-adu", 1, NULL, 3}, + {"start-params", 1, NULL, 4}, + {"push-res", 1, NULL, 5}, + {"res-push", 1, NULL, 5}, /* compat */ - {"no-scale", 0, &noscale, 1}, + {"no-scale", 0, &no_scale, 1}, {"no-polarisation", 0, &polarisation, 0}, {"no-polarization", 0, &polarisation, 0}, {"polarisation", 0, &polarisation, 1}, /* compat */ @@ -340,6 +554,20 @@ int main(int argc, char *argv[]) } break; + case 4 : + sparams_fn = strdup(optarg); + break; + + case 5 : + errno = 0; + push_res = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid value for --push-res.\n"); + return 1; + } + push_res = push_res*1e9; + break; + case 0 : break; @@ -392,12 +620,27 @@ int main(int argc, char *argv[]) } gsl_set_error_handler_off(); + rng = gsl_rng_alloc(gsl_rng_mt19937); /* Fill in what we know about the images so far */ n_images = 0; n_crystals = 0; images = NULL; crystals = NULL; + if ( sparams_fn != NULL ) { + char line[1024]; + sparams_fh = fopen(sparams_fn, "r"); + if ( sparams_fh == NULL ) { + ERROR("Failed to open '%s'\n", sparams_fn); + return 1; + } + fgets(line, 1024, sparams_fh); + STATUS("Reading initial scaling factors (G,B) from '%s'\n", + sparams_fn); + free(sparams_fn); + } else { + sparams_fh = NULL; + } do { @@ -460,10 +703,18 @@ int main(int argc, char *argv[]) cur); } + select_free_reflections(cr_refl, rng); + as = asymmetric_indices(cr_refl, sym); crystal_set_reflections(cr, as); + crystal_set_user_flag(cr, 0); reflist_free(cr_refl); + if ( set_initial_params(cr, sparams_fh) ) { + ERROR("Failed to set initial parameters\n"); + return 1; + } + n_crystals++; } @@ -475,6 +726,7 @@ int main(int argc, char *argv[]) } while ( 1 ); display_progress(n_images, n_crystals); fprintf(stderr, "\n"); + if ( sparams_fh != NULL ) fclose(sparams_fh); close_stream(st); @@ -484,102 +736,101 @@ int main(int argc, char *argv[]) for ( j=0; j<images[i].n_crystals; j++ ) { Crystal *cryst; - int n_gained = 0; - int n_lost = 0; - double mean_p_change = 0.0; cryst = images[i].crystals[j]; crystal_set_image(cryst, &images[i]); /* Now it's safe to do the following */ - update_partialities_2(cryst, pmodel, &n_gained, &n_lost, - &mean_p_change); - assert(n_gained == 0); /* That'd just be silly */ + update_partialities(cryst, pmodel); } } + /* Make a first pass at cutting out crap */ + STATUS("Checking patterns.\n"); + //early_rejection(crystals, n_crystals); + /* Make initial estimates */ - STATUS("Performing initial scaling.\n"); - if ( noscale ) STATUS("Scale factors fixed at 1.\n"); - full = scale_intensities(crystals, n_crystals, - nthreads, noscale, pmodel, min_measurements); + full = merge_intensities(crystals, n_crystals, nthreads, pmodel, + min_measurements, push_res); - srdata.crystals = crystals; - srdata.n = n_crystals; - srdata.full = full; - srdata.n_filtered = 0; - srdata.n_refined = 0; + check_rejection(crystals, n_crystals, full); - sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf", - infile, cmdline); - sr_iteration(sr, 0, &srdata); + show_duds(crystals, n_crystals); + + write_pgraph(full, crystals, n_crystals, 0); /* Iterate */ for ( i=0; i<n_iter; i++ ) { - int n_noscale = 0; - int n_noref = 0; - int n_solve = 0; - int n_lost = 0; - int n_dud = 0; - int j; - STATUS("Post refinement cycle %i of %i\n", i+1, n_iter); + double init_dev, init_free_dev; + double final_dev, final_free_dev; - srdata.n_filtered = 0; + STATUS("Refinement cycle %i of %i\n", i+1, n_iter); - /* Refine the geometry of all patterns to get the best fit */ + /* Refine all crystals to get the best fit */ refine_all(crystals, n_crystals, full, nthreads, pmodel, - &srdata); - - for ( j=0; j<n_crystals; j++ ) { - int flag; - flag = crystal_get_user_flag(crystals[j]); - if ( flag != 0 ) n_dud++; - if ( flag == 1 ) { - n_noscale++; - } else if ( flag == 2 ) { - n_noref++; - } else if ( flag == 3 ) { - n_solve++; - } else if ( flag == 4 ) { - n_lost++; - } - } + no_scale, &init_dev, &init_free_dev, + &final_dev, &final_free_dev); + + STATUS("Overall residual: initial = %e, final = %e\n", + init_dev, final_dev); + STATUS("Overall free residual: initial = %e, final = %e\n", + init_free_dev, final_free_dev); + + show_duds(crystals, n_crystals); + check_rejection(crystals, n_crystals, full); - if ( n_dud ) { - STATUS("%i crystals could not be refined this cycle.\n", - n_dud); - STATUS("%i scaling failed.\n", n_noscale); - STATUS("%i not enough reflections.\n", n_noref); - STATUS("%i solve failed.\n", n_solve); - STATUS("%i lost too many reflections.\n", n_lost); - } /* Re-estimate all the full intensities */ reflist_free(full); - full = scale_intensities(crystals, n_crystals, nthreads, - noscale, pmodel, min_measurements); - - srdata.full = full; + full = merge_intensities(crystals, n_crystals, nthreads, + pmodel, min_measurements, push_res); - sr_iteration(sr, i+1, &srdata); + write_pgraph(full, crystals, n_crystals, i+1); } - sr_finish(sr); - /* Output results */ write_reflist(outfile, full); + /* Output split results */ + char tmp[1024]; + RefList *split; + Crystal *crystals1[n_crystals]; + Crystal *crystals2[n_crystals]; + int n_crystals1 = 0; + int n_crystals2 = 0; + for ( i=0; i<n_crystals; i++ ) { + if ( i % 2 ) { + crystals1[n_crystals1] = crystals[i]; + n_crystals1++; + } else { + crystals2[n_crystals2] = crystals[i]; + n_crystals2++; + } + } + snprintf(tmp, 1024, "%s1", outfile); + split = merge_intensities(crystals1, n_crystals1, nthreads, + pmodel, min_measurements, push_res); + write_reflist(tmp, split); + reflist_free(split); + snprintf(tmp, 1024, "%s2", outfile); + split = merge_intensities(crystals2, n_crystals2, nthreads, + pmodel, min_measurements, push_res); + write_reflist(tmp, split); + reflist_free(split); + /* Dump parameters */ FILE *fh; fh = fopen("partialator.params", "w"); if ( fh == NULL ) { ERROR("Couldn't open partialator.params!\n"); } else { + fprintf(fh, " cr OSF relB div flag\n"); for ( i=0; i<n_crystals; i++ ) { - fprintf(fh, "%4i %5.2f %8.5e %s\n", i, + fprintf(fh, "%4i %10.5f %10.2f %8.5e %s\n", i, crystal_get_osf(crystals[i]), + crystal_get_Bfac(crystals[i])*1e20, crystal_get_image(crystals[i])->div, str_flags(crystals[i])); } @@ -587,6 +838,7 @@ int main(int argc, char *argv[]) } /* Clean up */ + gsl_rng_free(rng); for ( i=0; i<n_crystals; i++ ) { reflist_free(crystal_get_reflections(crystals[i])); crystal_free(crystals[i]); diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 4ee77c19..e37fd936 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -97,6 +97,7 @@ static void show_help(const char *s) " --beam-bandwidth Beam bandwidth as a fraction. Default 1%%.\n" " --photon-energy Photon energy in eV. Default 9000.\n" " --nphotons Number of photons per X-ray pulse. Default 1e12.\n" +" --beam-radius Radius of beam in metres (default 1e-6).\n" ); } @@ -661,7 +662,11 @@ int main(int argc, char *argv[]) STATUS("Simulation parameters:\n"); STATUS(" Photon energy: %.2f eV (wavelength %.5f A)\n", photon_energy, image.lambda*1e10); + STATUS(" Number of photons: %.0f (%.2f mJ)\n", nphotons, + eV_to_J(photon_energy)*1e3); STATUS(" Beam divergence: not simulated\n"); + STATUS(" Beam radius: %.2f microns\n", + beam_radius*1e6); STATUS(" Background: %.2f photons\n", background); switch ( spectrum_type ) { diff --git a/src/post-refinement.c b/src/post-refinement.c index 4dae5096..ccf010e2 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -54,7 +54,6 @@ /* Maximum number of iterations of NLSq to do for each image per macrocycle. */ #define MAX_CYCLES (10) - /* Returns dp(gauss)/dr at "r" */ static double gaussian_fraction_gradient(double r, double R) { @@ -175,31 +174,34 @@ static double volume_fraction(double rlow, double rhigh, double pr, } -/* Return the gradient of partiality wrt parameter 'k' given the current status - * of 'image'. */ -double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) +/* Return the gradient of "fx" wrt parameter 'k' given the current + * status of the crystal. */ +double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { - double azi; double glow, ghigh; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double xl, yl, zl; - double ds; - signed int hs, ks, ls; double rlow, rhigh, p; - double philow, phihigh, phi; - double khigh, klow; - double tl, cet, cez; struct image *image = crystal_get_image(cr); double R = crystal_get_profile_radius(cr); - double D, psph; + double gr; + signed int hi, ki, li; + double s; + get_indices(refl, &hi, &ki, &li); + s = resolution(crystal_get_cell(cr), hi, ki, li); get_partial(refl, &rlow, &rhigh, &p); - if ( k == REF_R ) { + if ( k == GPARAM_OSF ) { + return 1.0; + } + + if ( k == GPARAM_BFAC ) { + return -s*s; + } + + if ( k == GPARAM_R ) { double Rglow, Rghigh; + double D, psph; D = rlow - rhigh; psph = volume_fraction(rlow, rhigh, R, pmodel); @@ -207,7 +209,8 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) Rglow = volume_fraction_rgradient(rlow, R, pmodel); Rghigh = volume_fraction_rgradient(rhigh, R, pmodel); - return 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh); + gr = 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh); + return gr/p; } @@ -215,78 +218,23 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) glow = partiality_gradient(rlow, R, pmodel, rlow, rhigh); ghigh = partiality_gradient(rhigh, R, pmodel, rlow, rhigh); - get_symmetric_indices(refl, &hs, &ks, &ls); - ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); + if ( k == GPARAM_DIV ) { - cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - xl = hs*asx + ks*bsx + ls*csx; - yl = hs*asy + ks*bsy + ls*csy; - zl = hs*asz + ks*bsz + ls*csz; - - /* "low" gives the largest Ewald sphere (wavelength short => k large) - * "high" gives the smallest Ewald sphere (wavelength long => k small) - */ - klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); - khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0); - - tl = sqrt(xl*xl + yl*yl); - - cet = -sin(image->div/2.0) * klow; - cez = -cos(image->div/2.0) * klow; - philow = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0); - - cet = -sin(image->div/2.0) * khigh; - cez = -cos(image->div/2.0) * khigh; - phihigh = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0); - - /* Approximation: philow and phihigh are very similar */ - phi = (philow + phihigh) / 2.0; - - azi = atan2(yl, xl); - - /* For many gradients, just multiply the above number by the gradient - * of excitation error wrt whatever. */ - switch ( k ) { + double D, psph, ds; + signed int hs, ks, ls; - /* Cell parameters and orientation */ - case REF_ASX : - return - hs * sin(phi) * cos(azi) * (glow-ghigh); - - case REF_BSX : - return - ks * sin(phi) * cos(azi) * (glow-ghigh); - - case REF_CSX : - return - ls * sin(phi) * cos(azi) * (glow-ghigh); - - case REF_ASY : - return - hs * sin(phi) * sin(azi) * (glow-ghigh); - - case REF_BSY : - return - ks * sin(phi) * sin(azi) * (glow-ghigh); - - case REF_CSY : - return - ls * sin(phi) * sin(azi) * (glow-ghigh); - - case REF_ASZ : - return - hs * cos(phi) * (glow-ghigh); - - case REF_BSZ : - return - ks * cos(phi) * (glow-ghigh); - - case REF_CSZ : - return - ls * cos(phi) * (glow-ghigh); - - case REF_DIV : D = rlow - rhigh; psph = volume_fraction(rlow, rhigh, R, pmodel); - return (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D); + get_symmetric_indices(refl, &hs, &ks, &ls); + ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); + + gr = (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D); + return gr/p; } - ERROR("No gradient defined for parameter %i\n", k); - abort(); + gr = r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh); + return gr/p; } @@ -302,15 +250,15 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift) switch ( k ) { - case REF_ASX : asx += shift; break; - case REF_ASY : asy += shift; break; - case REF_ASZ : asz += shift; break; - case REF_BSX : bsx += shift; break; - case REF_BSY : bsy += shift; break; - case REF_BSZ : bsz += shift; break; - case REF_CSX : csx += shift; break; - case REF_CSY : csy += shift; break; - case REF_CSZ : csz += shift; break; + case GPARAM_ASX : asx += shift; break; + case GPARAM_ASY : asy += shift; break; + case GPARAM_ASZ : asz += shift; break; + case GPARAM_BSX : bsx += shift; break; + case GPARAM_BSY : bsy += shift; break; + case GPARAM_BSZ : bsz += shift; break; + case GPARAM_CSX : csx += shift; break; + case GPARAM_CSY : csy += shift; break; + case GPARAM_CSZ : csz += shift; break; } cell_set_reciprocal(cell, asx, asy, asz, @@ -325,32 +273,46 @@ static void apply_shift(Crystal *cr, int k, double shift) double t; struct image *image = crystal_get_image(cr); + if ( isnan(shift) ) { + ERROR("Refusing NaN shift for parameter %i\n", k); + ERROR("Image serial %i\n", image->serial); + return; + } + switch ( k ) { - case REF_DIV : - if ( isnan(shift) ) { - ERROR("NaN divergence shift\n"); - } else { - image->div += shift; - if ( image->div < 0.0 ) image->div = 0.0; - } + case GPARAM_DIV : + image->div += shift; + if ( image->div < 0.0 ) image->div = 0.0; break; - case REF_R : + case GPARAM_R : t = crystal_get_profile_radius(cr); t += shift; crystal_set_profile_radius(cr, t); break; - case REF_ASX : - case REF_ASY : - case REF_ASZ : - case REF_BSX : - case REF_BSY : - case REF_BSZ : - case REF_CSX : - case REF_CSY : - case REF_CSZ : + case GPARAM_BFAC : + t = crystal_get_Bfac(cr); + t += shift; + crystal_set_Bfac(cr, t); + break; + + case GPARAM_OSF : + t = crystal_get_osf(cr); + t += shift; + crystal_set_osf(cr, t); + break; + + case GPARAM_ASX : + case GPARAM_ASY : + case GPARAM_ASZ : + case GPARAM_BSX : + case GPARAM_BSY : + case GPARAM_BSZ : + case GPARAM_CSX : + case GPARAM_CSY : + case GPARAM_CSZ : apply_cell_shift(crystal_get_cell(cr), k, shift); break; @@ -362,143 +324,10 @@ static void apply_shift(Crystal *cr, int k, double shift) } -static int check_eigen(gsl_vector *e_val, int verbose) -{ - int i; - double vmax, vmin; - const int n = e_val->size; - const double max_condition = 1e6; - int n_filt = 0; - - if ( verbose ) STATUS("Eigenvalues:\n"); - vmin = +INFINITY; - vmax = 0.0; - for ( i=0; i<n; i++ ) { - double val = gsl_vector_get(e_val, i); - if ( verbose ) STATUS("%i: %e\n", i, val); - if ( val > vmax ) vmax = val; - if ( val < vmin ) vmin = val; - } - - for ( i=0; i<n; i++ ) { - double val = gsl_vector_get(e_val, i); - if ( val < vmax/max_condition ) { - gsl_vector_set(e_val, i, 0.0); - n_filt++; - } - } - - vmin = +INFINITY; - vmax = 0.0; - for ( i=0; i<n; i++ ) { - double val = gsl_vector_get(e_val, i); - if ( val == 0.0 ) continue; - if ( val > vmax ) vmax = val; - if ( val < vmin ) vmin = val; - } - if ( verbose ) { - STATUS("Condition number: %e / %e = %5.2f\n", - vmax, vmin, vmax/vmin); - STATUS("%i out of %i eigenvalues filtered.\n", n_filt, n); - } - - return n_filt; -} - - -static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt, - int verbose) -{ - gsl_matrix *s_vec; - gsl_vector *s_val; - int err, n; - gsl_vector *shifts; - gsl_vector *SB; - gsl_vector *SinvX; - gsl_matrix *S; /* rescaling matrix due to Bricogne */ - gsl_matrix *AS; - gsl_matrix *SAS; - int i; - gsl_matrix *SAS_copy; - - n = v->size; - if ( v->size != M->size1 ) return NULL; - if ( v->size != M->size2 ) return NULL; - - /* Calculate the rescaling matrix S */ - S = gsl_matrix_calloc(n, n); - for ( i=0; i<n; i++ ) { - double sii = pow(gsl_matrix_get(M, i, i), -0.5); - gsl_matrix_set(S, i, i, sii); - } - - /* Calculate the matrix SAS, which we will be (not) inverting */ - AS = gsl_matrix_calloc(n, n); - SAS = gsl_matrix_calloc(n, n); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, M, S, 0.0, AS); - gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, S, AS, 0.0, SAS); - gsl_matrix_free(AS); - - /* Calculate the vector SB, which is the RHS of the equation */ - SB = gsl_vector_calloc(n); - gsl_blas_dgemv(CblasNoTrans, 1.0, S, v, 0.0, SB); - - if ( verbose ) { - STATUS("The equation after rescaling:\n"); - show_matrix_eqn(SAS, SB); - } - - SAS_copy = gsl_matrix_alloc(SAS->size1, SAS->size2); - gsl_matrix_memcpy(SAS_copy, SAS); - - /* Do the SVD */ - s_val = gsl_vector_calloc(n); - s_vec = gsl_matrix_calloc(n, n); - err = gsl_linalg_SV_decomp_jacobi(SAS, s_vec, s_val); - if ( err ) { - if ( verbose ) ERROR("SVD failed: %s\n", gsl_strerror(err)); - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - gsl_matrix_free(SAS); - gsl_matrix_free(S); - return NULL; - } - /* "SAS" is now "U" */ - - /* Filter the eigenvalues */ - *n_filt = check_eigen(s_val, verbose); - - gsl_matrix_free(SAS_copy); - - /* Solve the equation SAS.SinvX = SB */ - SinvX = gsl_vector_calloc(n); - err = gsl_linalg_SV_solve(SAS, s_vec, s_val, SB, SinvX); - gsl_vector_free(SB); - gsl_matrix_free(SAS); - gsl_matrix_free(s_vec); - gsl_vector_free(s_val); - - if ( err ) { - ERROR("Matrix solution failed: %s\n", gsl_strerror(err)); - gsl_matrix_free(S); - gsl_vector_free(SinvX); - return NULL; - } - - /* Calculate S.SinvX to get X, the shifts */ - shifts = gsl_vector_calloc(n); - gsl_blas_dgemv(CblasNoTrans, 1.0, S, SinvX, 0.0, shifts); - - gsl_matrix_free(S); - gsl_vector_free(SinvX); - - return shifts; -} - - /* Perform one cycle of post refinement on 'image' against 'full' */ static double pr_iterate(Crystal *cr, const RefList *full, - PartialityModel pmodel, int *n_filtered) + PartialityModel pmodel, int no_scale, int *n_filtered, + int verbose) { gsl_matrix *M; gsl_vector *v; @@ -509,14 +338,42 @@ static double pr_iterate(Crystal *cr, const RefList *full, RefList *reflections; double max_shift; int nref = 0; - const int verbose = 0; + int num_params = 0; + enum gparam rv[32]; + double G, B; *n_filtered = 0; + /* If partiality model is anything other than "unity", refine all the + * geometrical parameters */ + if ( pmodel != PMODEL_UNITY ) { + rv[num_params++] = GPARAM_ASX; + rv[num_params++] = GPARAM_ASY; + rv[num_params++] = GPARAM_ASZ; + rv[num_params++] = GPARAM_BSX; + rv[num_params++] = GPARAM_BSY; + rv[num_params++] = GPARAM_BSZ; + rv[num_params++] = GPARAM_CSX; + rv[num_params++] = GPARAM_CSY; + rv[num_params++] = GPARAM_CSZ; + //rv[num_params++] = GPARAM_R; + } + + /* If we are scaling, refine scale factors (duh) */ + if ( !no_scale ) { + rv[num_params++] = GPARAM_OSF; + rv[num_params++] = GPARAM_BFAC; + } + + if ( num_params == 0 ) return 0.0; + reflections = crystal_get_reflections(cr); - M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS); - v = gsl_vector_calloc(NUM_PARAMS); + M = gsl_matrix_calloc(num_params, num_params); + v = gsl_vector_calloc(num_params); + + G = crystal_get_osf(cr); + B = crystal_get_Bfac(cr); /* Construct the equations, one per reflection in this image */ for ( refl = first_refl(reflections, &iter); @@ -524,54 +381,59 @@ static double pr_iterate(Crystal *cr, const RefList *full, refl = next_refl(refl, iter) ) { signed int ha, ka, la; - double I_full, delta_I; + double I_full, delta_I, esd; double w; double I_partial; int k; - double p, l; + double p, L, s; + double fx; Reflection *match; - double gradients[NUM_PARAMS]; + double gradients[num_params]; + + /* If reflection is free-flagged, don't use it here */ + if ( get_flag(refl) ) continue; /* Find the full version */ get_indices(refl, &ha, &ka, &la); match = find_refl(full, ha, ka, la); if ( match == NULL ) continue; - if ( (get_intensity(refl) < 3.0*get_esd_intensity(refl)) - || (get_partiality(refl) < MIN_PART_REFINE) - || (get_redundancy(match) < 2) ) continue; - + /* Merged intensitty */ I_full = get_intensity(match); - /* Actual measurement of this reflection from this pattern? */ - I_partial = get_intensity(refl) / crystal_get_osf(cr); + /* Actual measurement of this reflection from this pattern */ + I_partial = get_intensity(refl); + esd = get_esd_intensity(refl); + + if ( (get_partiality(refl) < MIN_PART_REFINE) + || (get_redundancy(match) < 2) + || (I_full <= 0) || (I_partial < 3*esd) ) continue; + p = get_partiality(refl); - l = get_lorentz(refl); + L = get_lorentz(refl); + s = resolution(crystal_get_cell(cr), ha, ka, la); /* Calculate the weight for this reflection */ - w = pow(get_esd_intensity(refl), 2.0); - w += l * p * I_full * pow(get_esd_intensity(match), 2.0); - w = pow(w, -1.0); + w = 1.0; /* Calculate all gradients for this reflection */ - for ( k=0; k<NUM_PARAMS; k++ ) { - gradients[k] = p_gradient(cr, k, refl, pmodel) * l; + for ( k=0; k<num_params; k++ ) { + gradients[k] = gradient(cr, rv[k], refl, pmodel); } - for ( k=0; k<NUM_PARAMS; k++ ) { + for ( k=0; k<num_params; k++ ) { int g; double v_c, v_curr; - for ( g=0; g<NUM_PARAMS; g++ ) { + for ( g=0; g<num_params; g++ ) { double M_c, M_curr; /* Matrix is symmetric */ if ( g > k ) continue; - M_c = gradients[g] * gradients[k]; - M_c *= w * pow(I_full, 2.0); + M_c = w * gradients[g] * gradients[k]; M_curr = gsl_matrix_get(M, k, g); gsl_matrix_set(M, k, g, M_curr + M_c); @@ -579,8 +441,9 @@ static double pr_iterate(Crystal *cr, const RefList *full, } - delta_I = I_partial - (l * p * I_full); - v_c = w * delta_I * I_full * gradients[k]; + fx = G + log(p) - log(L) - B*s*s + log(I_full); + delta_I = log(I_partial) - fx; + v_c = w * delta_I * gradients[k]; v_curr = gsl_vector_get(v, k); gsl_vector_set(v, k, v_curr + v_c); @@ -591,10 +454,10 @@ static double pr_iterate(Crystal *cr, const RefList *full, if ( verbose ) { STATUS("The original equation:\n"); show_matrix_eqn(M, v); + STATUS("%i reflections went into the equations.\n", nref); } - //STATUS("%i reflections went into the equations.\n", nref); - if ( nref == 0 ) { + if ( nref < num_params ) { crystal_set_user_flag(cr, 2); gsl_matrix_free(M); gsl_vector_free(v); @@ -605,10 +468,10 @@ static double pr_iterate(Crystal *cr, const RefList *full, shifts = solve_svd(v, M, n_filtered, verbose); if ( shifts != NULL ) { - for ( param=0; param<NUM_PARAMS; param++ ) { + for ( param=0; param<num_params; param++ ) { double shift = gsl_vector_get(shifts, param); - apply_shift(cr, param, shift); - //STATUS("Shift %i: %e\n", param, shift); + if ( verbose ) STATUS("Shift %i: %e\n", param, shift); + apply_shift(cr, rv[param], shift); if ( fabs(shift) > max_shift ) max_shift = fabs(shift); } @@ -624,107 +487,73 @@ static double pr_iterate(Crystal *cr, const RefList *full, } -static double guide_dev(Crystal *cr, const RefList *full) +static double residual(Crystal *cr, const RefList *full, int verbose, int free) { double dev = 0.0; - - /* For each reflection */ + double G, B; Reflection *refl; RefListIterator *iter; + FILE *fh = NULL; + + if ( verbose ) { + fh = fopen("residual.dat", "w"); + } + + G = crystal_get_osf(cr); + B = crystal_get_Bfac(cr); for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; - refl = next_refl(refl, iter) ) { - - double G, p; + refl = next_refl(refl, iter) ) + { + double p, L, s, w; signed int h, k, l; - Reflection *full_version; - double I_full, I_partial; + Reflection *match; + double esd, I_full, I_partial; + double fx, dc; - if ( (get_intensity(refl) < 3.0*get_esd_intensity(refl)) - || (get_partiality(refl) < MIN_PART_REFINE) ) continue; + if ( free && !get_flag(refl) ) continue; get_indices(refl, &h, &k, &l); - assert((h!=0) || (k!=0) || (l!=0)); - - full_version = find_refl(full, h, k, l); - if ( full_version == NULL ) continue; - /* Some reflections may have recently become scalable, but - * scale_intensities() might not yet have been called, so the - * full version may not have been calculated yet. */ + match = find_refl(full, h, k, l); + if ( match == NULL ) continue; - G = crystal_get_osf(cr); p = get_partiality(refl); + L = get_lorentz(refl); I_partial = get_intensity(refl); - I_full = get_intensity(full_version); - //STATUS("%3i %3i %3i %5.2f %5.2f %5.2f %5.2f %5.2f\n", - // h, k, l, G, p, I_partial, I_full, - // I_partial - p*G*I_full); - - dev += pow(I_partial - p*G*I_full, 2.0); - - } - - return dev; -} - - -struct param_backup -{ - UnitCell *cell; - double profile_radius; - double div; -}; - - -static struct param_backup backup_crystal(Crystal *cr) -{ - struct param_backup b; - struct image *image = crystal_get_image(cr); - - b.cell = cell_new_from_cell(crystal_get_cell(cr)); - b.profile_radius = crystal_get_profile_radius(cr); - b.div = image->div; - - return b; -} - - -static void revert_crystal(Crystal *cr, struct param_backup b) -{ - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - struct image *image = crystal_get_image(cr); + I_full = get_intensity(match); + esd = get_esd_intensity(refl); + s = resolution(crystal_get_cell(cr), h, k, l); - cell_get_reciprocal(b.cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); + if ( (get_partiality(refl) < MIN_PART_REFINE) + || (get_redundancy(match) < 2) + || (I_full <= 0) || (I_partial < 3*esd) ) continue; - cell_set_reciprocal(crystal_get_cell(cr), asx, asy, asz, - bsx, bsy, bsz, - csx, csy, csz); + fx = G + log(p) - log(L) - B*s*s + log(I_full); + dc = log(I_partial) - fx; + w = 1.0; + dev += w*dc*dc; - crystal_set_profile_radius(cr, b.profile_radius); - image->div = b.div; -} + if ( fh != NULL ) { + fprintf(fh, "%4i %4i %4i %e %.2f %e %f %f\n", + h, k, l, s, G, B, I_partial, I_full); + } + } + if ( fh != NULL ) fclose(fh); -static void free_backup_crystal(struct param_backup b) -{ - cell_free(b.cell); + return dev; } struct prdata pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel) + PartialityModel pmodel, int no_scale) { - double dev; int i; - struct param_backup backup; - const int verbose = 0; + int verbose = 0; struct prdata prdata; - double mean_p_change = 0.0; + int done = 0; + double old_dev; prdata.refined = 0; prdata.n_filtered = 0; @@ -732,15 +561,18 @@ struct prdata pr_refine(Crystal *cr, const RefList *full, /* Don't refine crystal if scaling was bad */ if ( crystal_get_user_flag(cr) != 0 ) return prdata; + old_dev = residual(cr, full, 0, 0); + prdata.initial_free_residual = residual(cr, full, 0, 1); + prdata.initial_residual = old_dev; + if ( verbose ) { - dev = guide_dev(cr, full); STATUS("\n"); /* Deal with progress bar */ - STATUS("Before iteration: dev = %10.5e\n", - dev); + STATUS("Initial G=%.2f, B=%e\n", + crystal_get_osf(cr), crystal_get_Bfac(cr)); + STATUS("Initial dev = %10.5e, free dev = %10.5e\n", + old_dev, residual(cr, full, 0, 1)); } - backup = backup_crystal(cr); - i = 0; do { @@ -748,44 +580,37 @@ struct prdata pr_refine(Crystal *cr, const RefList *full, double bsx, bsy, bsz; double csx, csy, csz; double dev; - int n_total; - int n_gained = 0; - int n_lost = 0; - n_total = num_reflections(crystal_get_reflections(cr)); cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - pr_iterate(cr, full, pmodel, &prdata.n_filtered); + pr_iterate(cr, full, pmodel, no_scale, &prdata.n_filtered, 0); - update_partialities_2(cr, pmodel, &n_gained, &n_lost, - &mean_p_change); + update_partialities(cr, pmodel); - if ( verbose ) { - dev = guide_dev(cr, full); - STATUS("PR Iteration %2i: mean p change = %10.2f" - " dev = %10.5e, %i gained, %i lost, %i total\n", - i+1, mean_p_change, dev, - n_gained, n_lost, n_total); - } + dev = residual(cr, full, 0, 0); + if ( fabs(dev - old_dev) < dev*0.01 ) done = 1; - if ( 3*n_lost > n_total ) { - revert_crystal(cr, backup); - update_partialities_2(cr, pmodel, &n_gained, &n_lost, - &mean_p_change); - crystal_set_user_flag(cr, 4); - break; + if ( verbose ) { + STATUS("Iter %2i: dev = %10.5e, free dev = %10.5e\n", + i+1, dev, residual(cr, full, 0, 1)); } i++; + old_dev = dev; - } while ( (mean_p_change > 0.01) && (i < MAX_CYCLES) ); - - free_backup_crystal(backup); + } while ( i < MAX_CYCLES && !done ); if ( crystal_get_user_flag(cr) == 0 ) { prdata.refined = 1; } + prdata.final_free_residual = residual(cr, full, 0, 1); + prdata.final_residual = old_dev; + + if ( verbose ) { + STATUS("Final G=%.2f, B=%e\n", crystal_get_osf(cr), + crystal_get_Bfac(cr)); + } return prdata; } diff --git a/src/post-refinement.h b/src/post-refinement.h index a9c79ed6..fd2d771b 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -43,36 +43,22 @@ #include "geometry.h" -/* Refineable parameters. - * Don't forget to add new things to backup_crystal() et al.! */ -enum { - REF_ASX, - REF_ASY, - REF_ASZ, - REF_BSX, - REF_BSY, - REF_BSZ, - REF_CSX, - REF_CSY, - REF_CSZ, - NUM_PARAMS, - REF_R, - REF_DIV, -}; - - struct prdata { int refined; int n_filtered; + double initial_residual; + double initial_free_residual; + double final_residual; + double final_free_residual; }; extern struct prdata pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel); + PartialityModel pmodel, int no_scale); /* Exported so it can be poked by tests/pr_p_gradient_check */ -extern double p_gradient(Crystal *cr, int k, Reflection *refl, - PartialityModel pmodel); +extern double gradient(Crystal *cr, int k, Reflection *refl, + PartialityModel pmodel); #endif /* POST_REFINEMENT_H */ diff --git a/src/predict-refine.c b/src/predict-refine.c new file mode 100644 index 00000000..2d845feb --- /dev/null +++ b/src/predict-refine.c @@ -0,0 +1,757 @@ +/* + * predict-refine.c + * + * Prediction refinement + * + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2015 Thomas White <taw@physics.org> + * + * 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 <stdlib.h> +#include <assert.h> +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> + +#include "image.h" +#include "geometry.h" +#include "cell-utils.h" + + +/* Maximum number of iterations of NLSq to do for each image per macrocycle. */ +#define MAX_CYCLES (10) + +/* Weighting of excitation error term (m^-1) compared to position term (m) */ +#define EXC_WEIGHT (4e-20) + +/* Parameters to refine */ +static const enum gparam rv[] = +{ + GPARAM_ASX, + GPARAM_ASY, + GPARAM_ASZ, + GPARAM_BSX, + GPARAM_BSY, + GPARAM_BSZ, + GPARAM_CSX, + GPARAM_CSY, + GPARAM_CSZ, + GPARAM_DETX, + GPARAM_DETY, +}; + +static const int num_params = 11; + +struct reflpeak { + Reflection *refl; + struct imagefeature *peak; + double Ih; /* normalised */ + struct panel *panel; /* panel the reflection appears on + * (we assume this never changes) */ +}; + + +static void twod_mapping(double fs, double ss, double *px, double *py, + struct panel *p) +{ + double xs, ys; + + xs = fs*p->fsx + ss*p->ssx; + ys = fs*p->fsy + ss*p->ssy; + + *px = (xs + p->cnx) / p->res; + *py = (ys + p->cny) / p->res; +} + + +static double r_dev(struct reflpeak *rp) +{ + /* Excitation error term */ + double rlow, rhigh, p; + get_partial(rp->refl, &rlow, &rhigh, &p); + return (rlow+rhigh)/2.0; +} + + +static double x_dev(struct reflpeak *rp, struct detector *det) +{ + /* Peak position term */ + double xpk, ypk, xh, yh; + double fsh, ssh; + twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel); + get_detector_pos(rp->refl, &fsh, &ssh); + twod_mapping(fsh, ssh, &xh, &yh, rp->panel); + return xh-xpk; +} + + +static double y_dev(struct reflpeak *rp, struct detector *det) +{ + /* Peak position term */ + double xpk, ypk, xh, yh; + double fsh, ssh; + twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel); + get_detector_pos(rp->refl, &fsh, &ssh); + twod_mapping(fsh, ssh, &xh, &yh, rp->panel); + return yh-ypk; +} + + +static void UNUSED write_pairs(const char *filename, struct reflpeak *rps, + int n, struct detector *det) +{ + int i; + FILE *fh; + + fh = fopen(filename, "w"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", filename); + return; + } + + for ( i=0; i<n; i++ ) { + + double write_fs, write_ss; + double fs, ss; + struct panel *p; + + fs = rps[i].peak->fs; + ss = rps[i].peak->ss; + + p = find_panel(det, fs, ss); + write_fs = fs - p->min_fs + p->orig_min_fs; + write_ss = ss - p->min_ss + p->orig_min_ss; + + fprintf(fh, "%7.2f %7.2f dev r,x,y: %9f %9f %9f %9f\n", + write_fs, write_ss, + r_dev(&rps[i])/1e9, fabs(r_dev(&rps[i])/1e9), + x_dev(&rps[i], det), + y_dev(&rps[i], det)); + + } + + fclose(fh); + + STATUS("Wrote %i pairs to %s\n", n, filename); +} + + +static int cmpd2(const void *av, const void *bv) +{ + struct reflpeak *a, *b; + + a = (struct reflpeak *)av; + b = (struct reflpeak *)bv; + + if ( fabs(r_dev(a)) < fabs(r_dev(b)) ) return -1; + return 1; +} + + +static int check_outlier_transition(struct reflpeak *rps, int n, + struct detector *det) +{ + int i; + + if ( n < 3 ) return n; + + qsort(rps, n, sizeof(struct reflpeak), cmpd2); + //write_pairs("pairs-before-outlier.lst", rps, n, det); + + for ( i=1; i<n-1; i++ ) { + + int j; + double grad = fabs(r_dev(&rps[i])) / i; + + for ( j=i+1; j<n; j++ ) { + if ( fabs(r_dev(&rps[j])) < 0.001e9+grad*j ) { + break; + } + } + if ( j == n ) { + //STATUS("Outlier transition found at position %i / %i\n", + // i, n); + return i; + } + } + + //STATUS("No outlier transition found.\n"); + return n; +} + + +/* Associate a Reflection with each peak in "image" which is close to Bragg. + * Reflections will be added to "reflist", which can be NULL if this is not + * needed. "rps" must be an array of sufficient size for all the peaks */ +static int pair_peaks(struct image *image, Crystal *cr, + RefList *reflist, struct reflpeak *rps) +{ + int i; + int n_acc = 0; + int n = 0; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + RefList *all_reflist; + + all_reflist = reflist_new(); + cell_get_cartesian(crystal_get_cell(cr), + &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + /* First, create a RefList containing the most likely indices for each + * peak, with no exclusion criteria */ + for ( i=0; i<image_feature_count(image->features); i++ ) { + + struct imagefeature *f; + double h, k, l, hd, kd, ld; + Reflection *refl; + struct panel *p; + + /* Assume all image "features" are genuine peaks */ + f = image_get_feature(image->features, i); + if ( f == NULL ) continue; + + /* Decimal and fractional Miller indices of nearest reciprocal + * lattice point */ + hd = f->rx * ax + f->ry * ay + f->rz * az; + kd = f->rx * bx + f->ry * by + f->rz * bz; + ld = f->rx * cx + f->ry * cy + f->rz * cz; + h = lrint(hd); + k = lrint(kd); + l = lrint(ld); + + refl = reflection_new(h, k, l); + if ( refl == NULL ) { + ERROR("Failed to create reflection\n"); + return 0; + } + + add_refl_to_list(refl, all_reflist); + set_symmetric_indices(refl, h, k, l); + + /* It doesn't matter if the actual predicted location + * doesn't fall on this panel. We're only interested + * in how far away it is from the peak location. + * The predicted position and excitation errors will be + * filled in by update_partialities(). */ + p = find_panel(image->det, f->fs, f->ss); + set_panel(refl, p); + + rps[n].refl = refl; + rps[n].peak = f; + rps[n].panel = p; + n++; + + } + + /* Get the excitation errors and detector positions for the candidate + * reflections */ + crystal_set_reflections(cr, all_reflist); + update_partialities(cr, PMODEL_SCSPHERE); + + /* Pass over the peaks again, keeping only the ones which look like + * good pairings */ + for ( i=0; i<n; i++ ) { + + double fs, ss, pd; + signed int h, k, l; + Reflection *refl = rps[i].refl; + + get_indices(refl, &h, &k, &l); + + /* Is the supposed reflection anywhere near the peak? */ + get_detector_pos(refl, &fs, &ss); + pd = pow(fs - rps[i].peak->fs, 2.0) + + pow(ss - rps[i].peak->ss, 2.0); + if ( pd > 10.0 * 10.0 ) continue; + + rps[n_acc] = rps[i]; + rps[n_acc].refl = reflection_new(h, k, l); + copy_data(rps[n_acc].refl, refl); + if ( reflist != NULL ) { + add_refl_to_list(rps[n_acc].refl, reflist); + } + n_acc++; + + } + reflist_free(all_reflist); + + /* Sort the pairings by excitation error and look for a transition + * between good pairings and outliers */ + n_acc = check_outlier_transition(rps, n_acc, image->det); + + return n_acc; +} + + +void refine_radius(Crystal *cr, struct image *image) +{ + int n, n_acc; + struct reflpeak *rps; + RefList *reflist; + + /* Maximum possible size */ + rps = malloc(image_feature_count(image->features) + * sizeof(struct reflpeak)); + if ( rps == NULL ) return; + + reflist = reflist_new(); + n_acc = pair_peaks(image, cr, reflist, rps); + if ( n_acc < 3 ) { + ERROR("Too few paired peaks (%i) to determine radius\n", n_acc); + free(rps); + return; + } + crystal_set_reflections(cr, reflist); + update_partialities(cr, PMODEL_SCSPHERE); + crystal_set_reflections(cr, NULL); + + qsort(rps, n_acc, sizeof(struct reflpeak), cmpd2); + n = (n_acc-1) - n_acc/50; + if ( n < 2 ) n = 2; /* n_acc is always >= 2 */ + crystal_set_profile_radius(cr, fabs(r_dev(&rps[n]))); + + reflist_free(reflist); + free(rps); +} + + +/* Returns d(xh-xpk)/dP, where P = any parameter */ +static double x_gradient(int param, struct reflpeak *rp, struct detector *det, + double lambda, UnitCell *cell) +{ + signed int h, k, l; + double xpk, ypk, xh, yh; + double fsh, ssh; + double x, z, wn; + double ax, ay, az, bx, by, bz, cx, cy, cz; + + twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel); + get_detector_pos(rp->refl, &fsh, &ssh); + twod_mapping(fsh, ssh, &xh, &yh, rp->panel); + get_indices(rp->refl, &h, &k, &l); + + wn = 1.0 / lambda; + + cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + x = h*ax + k*bx + l*cx; + z = h*az + k*bz + l*cz; + + switch ( param ) { + + case GPARAM_ASX : + return h * rp->panel->clen / (wn+z); + + case GPARAM_BSX : + return k * rp->panel->clen / (wn+z); + + case GPARAM_CSX : + return l * rp->panel->clen / (wn+z); + + case GPARAM_ASY : + return 0.0; + + case GPARAM_BSY : + return 0.0; + + case GPARAM_CSY : + return 0.0; + + case GPARAM_ASZ : + return -h * x * rp->panel->clen / (wn*wn + 2*wn*z + z*z); + + case GPARAM_BSZ : + return -k * x * rp->panel->clen / (wn*wn + 2*wn*z + z*z); + + case GPARAM_CSZ : + return -l * x * rp->panel->clen / (wn*wn + 2*wn*z + z*z); + + case GPARAM_DETX : + return -1; + + case GPARAM_DETY : + return 0; + + case GPARAM_CLEN : + return x / (wn+z); + + } + + ERROR("Positional gradient requested for parameter %i?\n", param); + abort(); +} + + +/* Returns d(yh-ypk)/dP, where P = any parameter */ +static double y_gradient(int param, struct reflpeak *rp, struct detector *det, + double lambda, UnitCell *cell) +{ + signed int h, k, l; + double xpk, ypk, xh, yh; + double fsh, ssh; + double y, z, wn; + double ax, ay, az, bx, by, bz, cx, cy, cz; + + twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel); + get_detector_pos(rp->refl, &fsh, &ssh); + twod_mapping(fsh, ssh, &xh, &yh, rp->panel); + get_indices(rp->refl, &h, &k, &l); + + wn = 1.0 / lambda; + + cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + y = h*ay + k*by + l*cy; + z = h*az + k*bz + l*cz; + + switch ( param ) { + + case GPARAM_ASX : + return 0.0; + + case GPARAM_BSX : + return 0.0; + + case GPARAM_CSX : + return 0.0; + + case GPARAM_ASY : + return h * rp->panel->clen / (wn+z); + + case GPARAM_BSY : + return k * rp->panel->clen / (wn+z); + + case GPARAM_CSY : + return l * rp->panel->clen / (wn+z); + + case GPARAM_ASZ : + return -h * y * rp->panel->clen / (wn*wn + 2*wn*z + z*z); + + case GPARAM_BSZ : + return -k * y * rp->panel->clen / (wn*wn + 2*wn*z + z*z); + + case GPARAM_CSZ : + return -l * y * rp->panel->clen / (wn*wn + 2*wn*z + z*z); + + case GPARAM_DETX : + return 0; + + case GPARAM_DETY : + return -1; + + case GPARAM_CLEN : + return y / (wn+z); + + } + + ERROR("Positional gradient requested for parameter %i?\n", param); + abort(); +} + + +static void update_detector(struct detector *det, double xoffs, double yoffs, + double coffs) +{ + int i; + + for ( i=0; i<det->n_panels; i++ ) { + struct panel *p = &det->panels[i]; + p->cnx += xoffs * p->res; + p->cny += yoffs * p->res; + p->clen += coffs; + } +} + + +static int iterate(struct reflpeak *rps, int n, UnitCell *cell, + struct image *image, + double *total_x, double *total_y, double *total_z) +{ + int i; + gsl_matrix *M; + gsl_vector *v; + gsl_vector *shifts; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + + /* Number of parameters to refine */ + M = gsl_matrix_calloc(num_params, num_params); + v = gsl_vector_calloc(num_params); + + for ( i=0; i<n; i++ ) { + + int k; + double gradients[num_params]; + double w; + + /* Excitation error terms */ + w = EXC_WEIGHT * rps[i].Ih; + + for ( k=0; k<num_params; k++ ) { + gradients[k] = r_gradient(cell, rv[k], rps[i].refl, + image); + } + + for ( k=0; k<num_params; k++ ) { + + int g; + double v_c, v_curr; + + for ( g=0; g<num_params; g++ ) { + + double M_c, M_curr; + + /* Matrix is symmetric */ + if ( g > k ) continue; + + M_c = w * gradients[g] * gradients[k]; + M_curr = gsl_matrix_get(M, k, g); + gsl_matrix_set(M, k, g, M_curr + M_c); + gsl_matrix_set(M, g, k, M_curr + M_c); + + } + + v_c = w * r_dev(&rps[i]); + v_c *= -gradients[k]; + v_curr = gsl_vector_get(v, k); + gsl_vector_set(v, k, v_curr + v_c); + + } + + /* Positional x terms */ + for ( k=0; k<num_params; k++ ) { + gradients[k] = x_gradient(rv[k], &rps[i], image->det, + image->lambda, cell); + } + + for ( k=0; k<num_params; k++ ) { + + int g; + double v_c, v_curr; + + for ( g=0; g<num_params; g++ ) { + + double M_c, M_curr; + + /* Matrix is symmetric */ + if ( g > k ) continue; + + M_c = gradients[g] * gradients[k]; + M_curr = gsl_matrix_get(M, k, g); + gsl_matrix_set(M, k, g, M_curr + M_c); + gsl_matrix_set(M, g, k, M_curr + M_c); + + } + + v_c = x_dev(&rps[i], image->det); + v_c *= -gradients[k]; + v_curr = gsl_vector_get(v, k); + gsl_vector_set(v, k, v_curr + v_c); + + } + + /* Positional y terms */ + for ( k=0; k<num_params; k++ ) { + gradients[k] = y_gradient(rv[k], &rps[i], image->det, + image->lambda, cell); + } + + for ( k=0; k<num_params; k++ ) { + + int g; + double v_c, v_curr; + + for ( g=0; g<num_params; g++ ) { + + double M_c, M_curr; + + /* Matrix is symmetric */ + if ( g > k ) continue; + + M_c = gradients[g] * gradients[k]; + M_curr = gsl_matrix_get(M, k, g); + gsl_matrix_set(M, k, g, M_curr + M_c); + gsl_matrix_set(M, g, k, M_curr + M_c); + + } + + v_c = y_dev(&rps[i], image->det); + v_c *= -gradients[k]; + v_curr = gsl_vector_get(v, k); + gsl_vector_set(v, k, v_curr + v_c); + + } + + } + + //show_matrix_eqn(M, v); + shifts = solve_svd(v, M, NULL, 0); + if ( shifts == NULL ) { + ERROR("Failed to solve equations.\n"); + gsl_matrix_free(M); + gsl_vector_free(v); + return 1; + } + + for ( i=0; i<num_params; i++ ) { + // STATUS("Shift %i = %e\n", i, gsl_vector_get(shifts, i)); + if ( isnan(gsl_vector_get(shifts, i)) ) { + gsl_vector_set(shifts, i, 0.0); + } + } + + /* Apply shifts */ + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + + /* Ensure the order here matches the order in rv[] */ + asx += gsl_vector_get(shifts, 0); + asy += gsl_vector_get(shifts, 1); + asz += gsl_vector_get(shifts, 2); + bsx += gsl_vector_get(shifts, 3); + bsy += gsl_vector_get(shifts, 4); + bsz += gsl_vector_get(shifts, 5); + csx += gsl_vector_get(shifts, 6); + csy += gsl_vector_get(shifts, 7); + csz += gsl_vector_get(shifts, 8); + update_detector(image->det, gsl_vector_get(shifts, 9), + gsl_vector_get(shifts, 10), + gsl_vector_get(shifts, 11)); + *total_x += gsl_vector_get(shifts, 9); + *total_y += gsl_vector_get(shifts, 10); + *total_z += gsl_vector_get(shifts, 11); + + cell_set_reciprocal(cell, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + + gsl_vector_free(shifts); + gsl_matrix_free(M); + gsl_vector_free(v); + + return 0; +} + + +static double UNUSED residual(struct reflpeak *rps, int n, struct detector *det) +{ + int i; + double res = 0.0; + double r; + + r = 0.0; + for ( i=0; i<n; i++ ) { + r += EXC_WEIGHT * rps[i].Ih * pow(r_dev(&rps[i]), 2.0); + } + printf("%e ", r); + res += r; + + r = 0.0; + for ( i=0; i<n; i++ ) { + r += pow(x_dev(&rps[i], det), 2.0); + } + printf("%e ", r); + res += r; + + r = 0.0; + for ( i=0; i<n; i++ ) { + r += pow(y_dev(&rps[i], det), 2.0); + } + printf("%e\n", r); + res += r; + + return res; +} + + +int refine_prediction(struct image *image, Crystal *cr) +{ + int n; + int i; + struct reflpeak *rps; + double max_I; + RefList *reflist; + double total_x = 0.0; + double total_y = 0.0; + double total_z = 0.0; + char tmp[1024]; + + rps = malloc(image_feature_count(image->features) + * sizeof(struct reflpeak)); + if ( rps == NULL ) return 1; + + reflist = reflist_new(); + n = pair_peaks(image, cr, reflist, rps); + if ( n < 10 ) { + ERROR("Too few paired peaks (%i) to refine orientation.\n", n); + free(rps); + reflist_free(reflist); + return 1; + } + crystal_set_reflections(cr, reflist); + + /* Normalise the intensities to max 1 */ + max_I = -INFINITY; + for ( i=0; i<n; i++ ) { + double cur_I = rps[i].peak->intensity; + if ( cur_I > max_I ) max_I = cur_I; + } + if ( max_I <= 0.0 ) { + ERROR("All peaks negative?\n"); + free(rps); + return 1; + } + for ( i=0; i<n; i++ ) { + rps[i].Ih = rps[i].peak->intensity / max_I; + } + + //STATUS("Initial residual = %e\n", residual(rps, n, image->det)); + + /* Refine */ + for ( i=0; i<MAX_CYCLES; i++ ) { + update_partialities(cr, PMODEL_SCSPHERE); + if ( iterate(rps, n, crystal_get_cell(cr), image, + &total_x, &total_y, &total_z) ) return 1; + //STATUS("Residual after %i = %e\n", i, + // residual(rps, n, image->det)); + } + //STATUS("Final residual = %e\n", residual(rps, n, image->det)); + + snprintf(tmp, 1024, "predict_refine/det_shift x = %.3f y = %.3f mm", + total_x*1e3, total_y*1e3); + crystal_add_notes(cr, tmp); + + crystal_set_reflections(cr, NULL); + reflist_free(reflist); + + n = pair_peaks(image, cr, NULL, rps); + free(rps); + if ( n < 10 ) { + ERROR("Too few paired peaks (%i) after refinement.\n", n); + return 1; + } + + return 0; +} diff --git a/src/predict-refine.h b/src/predict-refine.h new file mode 100644 index 00000000..c763d2ca --- /dev/null +++ b/src/predict-refine.h @@ -0,0 +1,45 @@ +/* + * predict-refine.h + * + * Prediction refinement + * + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2015 Thomas White <taw@physics.org> + * + * 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 PREDICT_REFINE_H +#define PREDICT_REFINE_H + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "crystal.h" + +struct image; + +extern int refine_prediction(struct image *image, Crystal *cr); +extern void refine_radius(Crystal *cr, struct image *image); + + +#endif /* PREDICT_REFINE_H */ diff --git a/src/process_hkl.c b/src/process_hkl.c index fbbf4ff0..6bddcafa 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -3,12 +3,12 @@ * * Assemble and process FEL Bragg intensities * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2014 Thomas White <taw@physics.org> + * 2009-2015 Thomas White <taw@physics.org> * 2011 Andrew Martin <andrew.martin@desy.de> * 2012 Lorenzo Galli <lorenzo.galli@desy.de> * 2014 Chunhong Yoon <chun.hong.yoon@desy.de> @@ -379,7 +379,6 @@ static void display_progress(int n_images, int n_crystals, int n_crystals_used) } - static int merge_all(Stream *st, RefList *model, RefList *reference, const SymOpList *sym, double **hist_vals, signed int hist_h, @@ -561,7 +560,8 @@ int main(int argc, char *argv[]) errno = 0; start_after = strtod(optarg, &rval); if ( *rval != '\0' ) { - ERROR("Invalid value for --start-after.\n"); + ERROR("Invalid value for --start-after (%s)\n", + optarg); return 1; } break; @@ -570,7 +570,8 @@ int main(int argc, char *argv[]) errno = 0; stop_after = strtod(optarg, &rval); if ( *rval != '\0' ) { - ERROR("Invalid value for --stop-after.\n"); + ERROR("Invalid value for --stop-after (%s)\n", + optarg); return 1; } break; diff --git a/src/process_image.c b/src/process_image.c index 6ccd0c9a..783ca464 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -3,11 +3,11 @@ * * The processing pipeline for one image * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * 2014 Valerio Mariani * * This file is part of CrystFEL. @@ -50,79 +50,35 @@ #include "reflist-utils.h" #include "process_image.h" #include "integration.h" +#include "predict-refine.h" -static int cmpd2(const void *av, const void *bv) +static void try_refine_autoR(struct image *image, Crystal *cr) { - double *ap, *bp; - double a, b; + double old_R, new_R; + char notes[1024]; - ap = (double *)av; - bp = (double *)bv; + refine_radius(cr, image); + old_R = crystal_get_profile_radius(cr); - a = ap[1]; - b = bp[1]; - - if ( fabs(a) < fabs(b) ) return -1; - return 1; -} - - -static void refine_radius(Crystal *cr) -{ - Reflection *refl; - RefListIterator *iter; - double vals[num_reflections(crystal_get_reflections(cr))*2]; - int n = 0; - int i; - double ti = 0.0; /* Total intensity */ - - for ( refl = first_refl(crystal_get_reflections(cr), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double rlow, rhigh, p; - int val = 0; - - if ( get_intensity(refl) > 9.0*get_esd_intensity(refl) ) { - val = 1; - } - - get_partial(refl, &rlow, &rhigh, &p); - - vals[(2*n)+0] = val; - vals[(2*n)+1] = fabs((rhigh+rlow)/2.0); - n++; - - } - - /* Sort in ascending order of absolute "deviation from Bragg" */ - qsort(vals, n, sizeof(double)*2, cmpd2); - - /* Calculate cumulative number of very strong reflections as a function - * of absolute deviation from Bragg */ - for ( i=0; i<n-1; i++ ) { - ti += vals[2*i]; - vals[2*i] = ti; - } - - if ( ti < 10 ) { - ERROR("WARNING: Not enough strong reflections (%.0f) to estimate " - "crystal parameters (trying anyway).\n", ti); + if ( refine_prediction(image, cr) ) { + crystal_set_user_flag(cr, 1); + return; } - /* Find the cutoff where we get 90% of the strong spots */ - for ( i=0; i<n-1; i++ ) { - if ( vals[2*i] > 0.90*ti ) break; - } + /* Estimate radius again with better geometry */ + refine_radius(cr, image); + new_R = crystal_get_profile_radius(cr); - crystal_set_profile_radius(cr, fabs(vals[2*i+1])); + snprintf(notes, 1024, "predict_refine/R old = %.5f new = %.5f nm^-1", + old_R/1e9, new_R/1e9); + crystal_add_notes(cr, notes); } void process_image(const struct index_args *iargs, struct pattern_args *pargs, Stream *st, int cookie, const char *tmpdir, int results_pipe, - int serial) + int serial, pthread_mutex_t *term_lock) { float *data_for_measurement; size_t data_size; @@ -133,6 +89,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, int r; int ret; char *rn; + int n_crystals_left; image.features = NULL; image.data = NULL; @@ -142,10 +99,11 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, image.filename = pargs->filename_p_e->filename; image.event = pargs->filename_p_e->ev; image.beam = iargs->beam; - image.det = iargs->det; + image.det = copy_geom(iargs->det); image.crystals = NULL; image.n_crystals = 0; image.serial = serial; + image.indexed_by = INDEXING_NONE; hdfile = hdfile_open(image.filename); if ( hdfile == NULL ) { @@ -177,10 +135,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, switch ( iargs->peaks ) { case PEAK_HDF5: - /* Get peaks from HDF5 */ - - if ( get_peaks(&image, hdfile, iargs->hdf5_peak_path, - iargs->cxi_hdf5_peaks, pargs->filename_p_e) ) { + if ( get_peaks(&image, hdfile, iargs->hdf5_peak_path) ) { ERROR("Failed to get peaks from HDF5 file.\n"); } if ( !iargs->no_revalidate ) { @@ -191,6 +146,19 @@ 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) ) { + ERROR("Failed to get peaks from CXI file.\n"); + } + if ( !iargs->no_revalidate ) { + validate_peaks(&image, iargs->min_snr, + iargs->pk_inn, iargs->pk_mid, + iargs->pk_out, iargs->use_saturated, + iargs->check_hdf5_snr); + } + break; + case PEAK_ZAEF: search_peaks(&image, iargs->threshold, iargs->min_gradient, iargs->min_snr, @@ -226,36 +194,73 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, } free(rn); - pargs->n_crystals = image.n_crystals; for ( i=0; i<image.n_crystals; i++ ) { crystal_set_image(image.crystals[i], &image); + crystal_set_user_flag(image.crystals[i], 0); } - /* Default parameters */ - image.div = 0.0; - image.bw = 0.00000001; - for ( i=0; i<image.n_crystals; i++ ) { - crystal_set_profile_radius(image.crystals[i], 0.01e9); - crystal_set_mosaicity(image.crystals[i], 0.0); /* radians */ + /* Set beam/crystal parameters */ + if ( iargs->fix_divergence >= 0.0 ) { + image.div = iargs->fix_divergence; + } else { + image.div = 0.0; + } + if ( iargs->fix_bandwidth >= 0.0 ) { + image.bw = iargs->fix_bandwidth; + } else { + image.bw = 0.00000001; + } + if ( iargs->fix_profile_r >= 0.0 ) { + for ( i=0; i<image.n_crystals; i++ ) { + crystal_set_profile_radius(image.crystals[i], + iargs->fix_profile_r); + crystal_set_mosaicity(image.crystals[i], 0.0); + } + } else { + for ( i=0; i<image.n_crystals; i++ ) { + crystal_set_profile_radius(image.crystals[i], 0.02e9); + crystal_set_mosaicity(image.crystals[i], 0.0); + } } - /* Integrate all the crystals at once - need all the crystals so that - * overlaps can be detected. */ - integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE, iargs->push_res, - iargs->ir_inn, iargs->ir_mid, iargs->ir_out, - INTDIAG_NONE, 0, 0, 0, results_pipe); + if ( iargs->fix_profile_r < 0.0 ) { + + for ( i=0; i<image.n_crystals; i++ ) { + if ( iargs->predict_refine ) { + try_refine_autoR(&image, image.crystals[i]); + } else { + refine_radius(image.crystals[i], &image); + } + } + + } else { + for ( i=0; i<image.n_crystals; i++ ) { + if ( iargs->predict_refine ) { + refine_prediction(&image, image.crystals[i]); + } + } + + } + + /* If there are no crystals left, set the indexing flag back to zero */ + n_crystals_left = 0; for ( i=0; i<image.n_crystals; i++ ) { - refine_radius(image.crystals[i]); - reflist_free(crystal_get_reflections(image.crystals[i])); + if ( crystal_get_user_flag(image.crystals[i]) == 0 ) { + n_crystals_left++; + } + } + if ( n_crystals_left == 0 ) { + image.indexed_by = INDEXING_NONE; } + /* Integrate! */ integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE, - iargs->push_res, - iargs->ir_inn, iargs->ir_mid, iargs->ir_out, - iargs->int_diag, iargs->int_diag_h, - iargs->int_diag_k, iargs->int_diag_l, - results_pipe); + iargs->push_res, + iargs->ir_inn, iargs->ir_mid, iargs->ir_out, + iargs->int_diag, iargs->int_diag_h, + iargs->int_diag_k, iargs->int_diag_l, + term_lock); ret = write_chunk(st, &image, hdfile, iargs->stream_peaks, iargs->stream_refls, @@ -269,8 +274,17 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, n += crystal_get_num_implausible_reflections(image.crystals[i]); } if ( n > 0 ) { - STATUS("WARNING: %i implausibly negative reflection%s in %s.\n", - n, n>1?"s":"", image.filename); + STATUS("WARNING: %i implausibly negative reflection%s in %s " + "%s\n", n, n>1?"s":"", image.filename, + get_event_string(image.event)); + } + + /* Count crystals which are still good */ + pargs->n_crystals = 0; + for ( i=0; i<image.n_crystals; i++ ) { + if ( crystal_get_user_flag(image.crystals[i]) == 0 ) { + pargs->n_crystals++; + } } for ( i=0; i<image.n_crystals; i++ ) { @@ -290,5 +304,6 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, free(image.data); if ( image.flags != NULL ) free(image.flags); image_feature_list_free(image.features); + free_detector_geometry(image.det); hdfile_close(hdfile); } diff --git a/src/process_image.h b/src/process_image.h index a0bd6a83..d982c4f0 100644 --- a/src/process_image.h +++ b/src/process_image.h @@ -3,11 +3,11 @@ * * The processing pipeline for one image * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2014 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * 2014 Valerio Mariani * * This file is part of CrystFEL. @@ -41,6 +41,7 @@ enum { PEAK_ZAEF, PEAK_HDF5, + PEAK_CXI, }; @@ -63,7 +64,6 @@ struct index_args float tols[4]; struct beam_params *beam; char *hdf5_peak_path; - int cxi_hdf5_peaks; float pk_inn; float pk_mid; float pk_out; @@ -83,6 +83,10 @@ struct index_args signed int int_diag_l; float push_res; float highres; + float fix_profile_r; + float fix_bandwidth; + float fix_divergence; + int predict_refine; }; @@ -100,7 +104,7 @@ struct pattern_args extern void process_image(const struct index_args *iargs, struct pattern_args *pargs, Stream *st, int cookie, const char *tmpdir, int results_pipe, - int serial); + int serial, pthread_mutex_t *term_lock); #endif /* PROCESS_IMAGEs_H */ diff --git a/src/rejection.c b/src/rejection.c new file mode 100644 index 00000000..a565fec3 --- /dev/null +++ b/src/rejection.c @@ -0,0 +1,191 @@ +/* + * rejection.c + * + * Crystal rejection for scaling + * + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2015 Thomas White <taw@physics.org> + * + * 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 <stdlib.h> +#include <assert.h> +#include <gsl/gsl_statistics.h> + +#include "crystal.h" +#include "reflist.h" +#include "rejection.h" +#include "cell-utils.h" + + +static double mean_intensity(RefList *list) +{ + Reflection *refl; + RefListIterator *iter; + double total = 0.0; + int n = 0; + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + total += get_intensity(refl); + n++; + } + + return total/n; +} + + +/* Reject really obvious outliers */ +void early_rejection(Crystal **crystals, int n) +{ + int i; + double m = 0.0; + double mean_m; + FILE *fh = fopen("reject.dat", "w"); + int n_flag = 0; + + for ( i=0; i<n; i++ ) { + double u; + RefList *list = crystal_get_reflections(crystals[i]); + u = mean_intensity(list); + m += u; + fprintf(fh, "%i %f\n", i, u); + } + mean_m = m/n; + for ( i=0; i<n; i++ ) { + double u; + RefList *list = crystal_get_reflections(crystals[i]); + u = mean_intensity(list); + if ( u/mean_m < 0.2 ) { + crystal_set_user_flag(crystals[i], 5); + n_flag++; + } + } + fclose(fh); + + STATUS("Mean intensity/peak = %f ADU\n", m/n); + STATUS("%i crystals flagged\n", n_flag); +} + + +static void check_cc(Crystal *cr, RefList *full) +{ + RefList *list = crystal_get_reflections(cr); + Reflection *refl; + RefListIterator *iter; + double G = crystal_get_osf(cr); + double B = crystal_get_Bfac(cr); + UnitCell *cell = crystal_get_cell(cr); + double *vec1, *vec2; + int n, i; + double cc; + + n = num_reflections(list); + vec1 = malloc(n*sizeof(double)); + vec2 = malloc(n*sizeof(double)); + if ( (vec1 == NULL) || (vec2 == NULL) ) { + ERROR("Not enough memory to check CCs\n"); + return; + } + + i = 0; + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + double pobs, pcalc; + double res, corr, Ipart; + Reflection *match; + + if ( !get_flag(refl) ) continue; /* Not free-flagged */ + + get_indices(refl, &h, &k, &l); + res = resolution(cell, h, k, l); + if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue; + + match = find_refl(full, h, k, l); + if ( match == NULL ) continue; + + /* Calculated partiality */ + pcalc = get_partiality(refl); + + /* Observed partiality */ + corr = exp(-G) * exp(B*res*res) * get_lorentz(refl); + Ipart = get_intensity(refl) * corr; + pobs = Ipart / get_intensity(match); + + vec1[i] = pobs; + vec2[i] = pcalc; + i++; + } + + cc = gsl_stats_correlation(vec1, 1, vec2, 1, i); + //printf("%f\n", cc); + if ( cc < 0.5 ) crystal_set_user_flag(cr, 6); + + free(vec1); + free(vec2); +} + + +static void check_ccs(Crystal **crystals, int n_crystals, RefList *full) +{ + int i; + + for ( i=0; i<n_crystals; i++ ) { + check_cc(crystals[i], full); + } +} + + +void check_rejection(Crystal **crystals, int n, RefList *full) +{ + int i; + int n_acc = 0; + + /* Check according to CCs FIXME: Disabled */ + //if ( full != NULL ) check_ccs(crystals, n, full); + + for ( i=0; i<n; i++ ) { + + /* Reject if B factor modulus is very large */ + if ( fabs(crystal_get_Bfac(crystals[i])) > 1e-18 ) { + crystal_set_user_flag(crystals[i], 1); + } + + if ( crystal_get_user_flag(crystals[i]) == 0 ) n_acc++; + + } + + if ( n_acc < 2 ) { + ERROR("Not enough crystals left to proceed (%i). Sorry.\n", + n_acc); + exit(1); + } +} diff --git a/src/rejection.h b/src/rejection.h new file mode 100644 index 00000000..ec529941 --- /dev/null +++ b/src/rejection.h @@ -0,0 +1,43 @@ +/* + * rejection.h + * + * Crystal rejection for scaling + * + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2015 Thomas White <taw@physics.org> + * + * 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 REJECTION_H +#define REJECTION_H + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include "crystal.h" + +extern void early_rejection(Crystal **crystals, int n); +extern void check_rejection(Crystal **crystals, int n, RefList *full); + +#endif /* REJECTION_H */ diff --git a/src/scaling-report.c b/src/scaling-report.c deleted file mode 100644 index ca4c5cbc..00000000 --- a/src/scaling-report.c +++ /dev/null @@ -1,898 +0,0 @@ -/* - * scaling-report.c - * - * Write a nice PDF of scaling parameters - * - * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2010-2013 Thomas White <taw@physics.org> - * - * 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 <cairo.h> -#include <cairo-pdf.h> -#include <pango/pangocairo.h> -#include <math.h> - -#include "image.h" -#include "scaling-report.h" - - -#define PAGE_WIDTH (842.0) - -enum justification -{ - J_CENTER, - J_LEFT, - J_RIGHT, -}; - - -struct _srcontext -{ - cairo_surface_t *surf; - cairo_t *cr; - double w; - double h; - - /* Most sampled reflections */ - signed int ms_h[9]; - signed int ms_k[9]; - signed int ms_l[9]; - -}; - - -static void show_text(cairo_t *cr, const char *text, double y, - enum justification j, char *font) -{ - PangoLayout *layout; - PangoFontDescription *fontdesc; - int width, height; - PangoAlignment just; - - if ( font == NULL ) font = "Sans 10"; - - layout = pango_cairo_create_layout(cr); - pango_layout_set_ellipsize(layout, PANGO_ELLIPSIZE_NONE); - pango_layout_set_width(layout, PANGO_SCALE*(PAGE_WIDTH-20.0)); - - switch ( j ) - { - case J_CENTER : just = PANGO_ALIGN_CENTER; break; - case J_LEFT : just = PANGO_ALIGN_LEFT; break; - case J_RIGHT : just = PANGO_ALIGN_RIGHT; break; - default: just = PANGO_ALIGN_LEFT; break; - } - - pango_layout_set_alignment(layout, just); - pango_layout_set_wrap(layout, PANGO_WRAP_CHAR); - pango_layout_set_spacing(layout, 4.0*PANGO_SCALE); - - pango_layout_set_text(layout, text, -1); - - fontdesc = pango_font_description_from_string(font); - pango_layout_set_font_description(layout, fontdesc); - - pango_cairo_update_layout(cr, layout); - pango_layout_get_size(layout, &width, &height); - - cairo_move_to(cr, 10.0, y); - cairo_set_source_rgb(cr, 0.0, 0.0, 0.0); - pango_cairo_show_layout(cr, layout); -} - - -static void show_text_simple(cairo_t *cr, const char *text, double x, double y, - char *font, double rot, enum justification j) -{ - PangoLayout *layout; - PangoFontDescription *fontdesc; - int width, height; - - cairo_save(cr); - - if ( font == NULL ) font = "Sans 10"; - - layout = pango_cairo_create_layout(cr); - pango_layout_set_ellipsize(layout, PANGO_ELLIPSIZE_NONE); - pango_layout_set_alignment(layout, PANGO_ALIGN_LEFT); - - pango_layout_set_text(layout, text, -1); - - fontdesc = pango_font_description_from_string(font); - pango_layout_set_font_description(layout, fontdesc); - - pango_cairo_update_layout(cr, layout); - pango_layout_get_size(layout, &width, &height); - - cairo_new_path(cr); - cairo_translate(cr, x, y); - cairo_rotate(cr, rot); - if ( j == J_CENTER ) { - cairo_translate(cr, -(width/2.0)/PANGO_SCALE, - -(height/2.0)/PANGO_SCALE); - } else if ( j == J_RIGHT ) { - cairo_translate(cr, -width/PANGO_SCALE, - -(height/2.0)/PANGO_SCALE); - } else { - cairo_translate(cr, 0.0, -(height/2.0)/PANGO_SCALE); - } - cairo_set_source_rgb(cr, 0.0, 0.0, 0.0); - pango_cairo_show_layout(cr, layout); - - cairo_restore(cr); -} - - -static void plot_point(cairo_t *cr, double g_width, double g_height, - double pcalc, double pobs) -{ - int bad = 0; - - if ( pobs > 1.0 ) { - pobs = 1.01; - bad = 1; - - } - if ( pcalc > 1.0 ) { - pcalc = 1.01; - bad = 1; - } - if ( pobs < 0.0 ) { - pobs = -0.01; - bad = 1; - } - if ( pcalc < 0.0 ) { - pobs = -0.01; - bad = 1; - } - - if ( bad ) { - cairo_set_source_rgb(cr, 1.0, 0.0, 0.0); - } - - cairo_arc(cr, g_width*pcalc, g_height*(1.0-pobs), - 1.0, 0.0, 2.0*M_PI); - cairo_fill(cr); - - if ( bad ) { - cairo_set_source_rgb(cr, 0.0, 0.7, 0.0); - } -} - - -static void partiality_graph(cairo_t *cr, Crystal **crystals, int n, - RefList *full) -{ - const double g_width = 200.0; - const double g_height = 200.0; - int i; - const int nbins = 25; - double t_num[nbins]; - double t_den[nbins]; - double prob; - double pcalcmin[nbins]; - double pcalcmax[nbins]; - int num_nondud; - gsl_rng *rng; - - show_text_simple(cr, "Observed partiality", -20.0, g_height/2.0, - NULL, -M_PI_2, J_CENTER); - show_text_simple(cr, "Calculated partiality", - g_width/2.0,g_height+20.0, NULL, 0.0, J_CENTER); - - show_text_simple(cr, "0.0", -20.0, g_height, NULL, 0.0, J_CENTER); - show_text_simple(cr, "1.0", -20.0, 0.0, NULL, 0.0, J_CENTER); - show_text_simple(cr, "0.0", 0.0, g_height+10.0, NULL, - -M_PI/3.0, J_RIGHT); - show_text_simple(cr, "1.0", g_width, g_height+10.0, NULL, - -M_PI/3.0, J_RIGHT); - - for ( i=0; i<nbins; i++ ) { - t_num[i] = 0.0; - t_den[i] = 0.0; - pcalcmin[i] = (double)i/nbins; - pcalcmax[i] = (double)(i+1)/nbins; - } - pcalcmax[nbins-1] += 0.001; /* Make sure it include pcalc = 1 */ - - num_nondud = 0; - for ( i=0; i<n; i++ ) { - if ( crystal_get_user_flag(crystals[i]) ) continue; - num_nondud++; - } - - /* The reflections chosen for the graph will be the same every time - * (given the same sequence of input reflections, scalabilities etc) */ - rng = gsl_rng_alloc(gsl_rng_mt19937); - - cairo_set_source_rgb(cr, 0.0, 0.7, 0.0); - prob = 1.0 / num_nondud; - for ( i=0; i<n; i++ ) { - - Reflection *refl; - RefListIterator *iter; - Crystal *cryst; - - cryst = crystals[i]; - - if ( crystal_get_user_flag(cryst) ) continue; - - for ( refl = first_refl(crystal_get_reflections(cryst), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double Ipart, Ifull, pobs, pcalc; - signed int h, k, l; - Reflection *f; - int bin; - - get_indices(refl, &h, &k, &l); - f = find_refl(full, h, k, l); - if ( f == NULL ) continue; - if ( get_redundancy(f) < 2 ) continue; - - Ipart = crystal_get_osf(cryst) * get_intensity(refl); - Ifull = get_intensity(f); - - pobs = Ipart/Ifull; - pcalc = get_lorentz(refl) * get_partiality(refl); - - //STATUS("%4i %4i %4i : %9.6f %9.6f %e %e %e\n", - // h, k, l, pobs, pcalc, - // Ipart, Ifull, images[i].osf); - - for ( bin=0; bin<nbins; bin++ ) { - if ( (pcalc >= pcalcmin[bin]) - && (pcalc < pcalcmax[bin]) ) - { - double esd_pobs, esd_Ip, esd_If; - esd_Ip = get_esd_intensity(refl); - esd_If = get_esd_intensity(f); - esd_If *= crystal_get_osf(cryst); - esd_pobs = pow(esd_Ip/Ipart, 2.0); - esd_pobs += pow(esd_If/Ifull, 2.0); - esd_pobs = sqrt(esd_pobs); - t_num[bin] += pobs / esd_pobs; - t_den[bin] += 1.0 / esd_pobs; - } - } - - bin = nbins * pcalc; - - if ( random_flat(rng, 1.0) < prob ) { - plot_point(cr, g_width, g_height, pcalc, pobs); - } - } - - } - - gsl_rng_free(rng); - - cairo_new_path(cr); - cairo_rectangle(cr, 0.0, 0.0, g_width, g_height); - cairo_clip(cr); - - cairo_new_path(cr); - cairo_move_to(cr, 0.0, g_height); - for ( i=0; i<nbins; i++ ) { - - double pos = pcalcmin[i] + (pcalcmax[i] - pcalcmin[i])/2.0; - - if ( t_den[i] == 0.0 ) continue; - cairo_line_to(cr, g_width*pos, - g_height - g_height*(t_num[i]/t_den[i])); - - } - cairo_set_source_rgb(cr, 0.0, 0.0, 0.0); - cairo_stroke(cr); - - cairo_reset_clip(cr); - - cairo_new_path(cr); - cairo_rectangle(cr, 0.0, 0.0, g_width, g_height); - cairo_set_source_rgb(cr, 0.0, 0.0, 0.0); - cairo_set_line_width(cr, 1.5); - cairo_stroke(cr); -} - - -static void partiality_histogram(cairo_t *cr, Crystal **crystals, int n, - RefList *full, int calc, int backwards) -{ - int f_max; - int i, b; - const int nbins = 100; - int counts[nbins]; - const double g_width = 200.0; - const double g_height = 120.0; - char tmp[32]; - double text_rot, axis_pos; - - if ( backwards ) { - text_rot = M_PI_2; - axis_pos = 215.0; - } else { - text_rot = -M_PI_2; - axis_pos = -15.0; - } - - show_text_simple(cr, "Frequency", axis_pos, g_height/2.0, - NULL, text_rot, J_CENTER); - - for ( b=0; b<nbins; b++ ) { - counts[b] = 0; - } - - for ( i=0; i<n; i++ ) { - - Reflection *refl; - RefListIterator *iter; - Crystal *cryst = crystals[i]; - - if ( crystal_get_user_flag(cryst) ) continue; - - for ( refl = first_refl(crystal_get_reflections(cryst), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double Ipart, Ifull, pobs, pcalc; - signed int h, k, l; - Reflection *f; - - get_indices(refl, &h, &k, &l); - f = find_refl(full, h, k, l); - if ( f == NULL ) continue; - - Ipart = get_intensity(refl); - Ifull = get_intensity(f); - - pobs = (Ipart * crystal_get_osf(cryst)) - / (Ifull * get_lorentz(refl)); - pcalc = get_partiality(refl); - - if ( calc ) { - b = pcalc*nbins; - } else { - b = pobs*nbins; - } - if ( (b>=0) && (b<nbins) ) counts[b]++; - } - } - - f_max = 0; - for ( b=0; b<nbins; b++ ) { - if ( counts[b] > f_max ) f_max = counts[b]; - } - f_max = (f_max/10)*10 + 10; - - if ( !backwards ) { - show_text_simple(cr, "0", axis_pos, g_height, - NULL, 0.0, J_RIGHT); - } else { - show_text_simple(cr, "0", axis_pos, 0.0, - NULL, M_PI, J_RIGHT); - } - snprintf(tmp, 31, "%i", f_max); - if ( !backwards ) { - show_text_simple(cr, tmp, axis_pos, 0.0, - NULL, 0.0, J_RIGHT); - } else { - show_text_simple(cr, tmp, axis_pos, g_height, - NULL, M_PI, J_RIGHT); - } - - for ( b=0; b<nbins; b++ ) { - - double bar_height; - - bar_height = ((double)counts[b]/f_max)*g_height; - - cairo_new_path(cr); - if ( !backwards ) { - cairo_rectangle(cr, (g_width/nbins)*b, g_height, - g_width/nbins, -bar_height); - } else { - cairo_rectangle(cr, (g_width/nbins)*b, 0.0, - g_width/nbins, bar_height); - } - cairo_set_source_rgb(cr, 0.0, 0.0, 1.0); - cairo_set_line_width(cr, 1.0); - cairo_stroke(cr); - - } - - cairo_new_path(cr); - cairo_rectangle(cr, 0.0, 0.0, g_width, g_height); - cairo_set_source_rgb(cr, 0.0, 0.0, 0.0); - cairo_set_line_width(cr, 1.5); - cairo_stroke(cr); -} - - -static void scale_factor_histogram(cairo_t *cr, Crystal **crystals, int n, - const char *title) -{ - int f_max; - int i, b; - const int nbins = 100; - double osf_max, osf_inc; - double osf_low[nbins]; - double osf_high[nbins]; - int counts[nbins]; - const double g_width = 320.0; - const double g_height = 180.0; - char tmp[32]; - int n_zero, n_half; - - show_text_simple(cr, title, g_width/2.0, -18.0, - "Sans Bold 10", 0.0, J_CENTER); - - show_text_simple(cr, "Frequency", -15.0, g_height/2.0, - NULL, -M_PI_2, J_CENTER); - show_text_simple(cr, "Scale factor", g_width/2.0, g_height+12.0, - NULL, 0.0, J_CENTER); - - osf_max = 0.0; - for ( i=0; i<n; i++ ) { - double osf = crystal_get_osf(crystals[i]); - if ( crystal_get_user_flag(crystals[i]) ) continue; - if ( osf > osf_max ) osf_max = osf; - } - osf_max = ceil(osf_max+osf_max/10000.0); - if ( osf_max > 1000.0 ) { - ERROR("Silly scale factor detected. Using 100.0 instead.\n"); - osf_max = 100.0; - } - - do { - - osf_inc = osf_max / nbins; - - for ( b=0; b<nbins; b++ ) { - osf_low[b] = b*osf_inc; - osf_high[b] = (b+1)*osf_inc; - counts[b] = 0; - } - - for ( i=0; i<n; i++ ) { - - double osf = crystal_get_osf(crystals[i]); - - if ( crystal_get_user_flag(crystals[i]) ) continue; - - for ( b=0; b<nbins; b++ ) { - if ( (osf >= osf_low[b]) - && (osf < osf_high[b]) ) - { - counts[b]++; - break; - } - } - } - - n_zero = 0; - n_half = 0; - if ( osf_max > 10.0 ) { - - /* Count the number of bins with no counts, subtract 1 - * from the maximum value until this isn't the case */ - for ( b=0; b<nbins; b++ ) { - if ( counts[b] == 0 ) n_zero++; - n_half++; - } - n_half /= 2; - - if ( n_zero > n_half ) osf_max -= 1.0; - - } - - } while ( n_zero > n_half ); - - f_max = 0; - for ( b=0; b<nbins; b++ ) { - if ( counts[b] > f_max ) f_max = counts[b]; - } - f_max = (f_max/10)*10 + 10; - - show_text_simple(cr, "0", -10.0, g_height, NULL, 0.0, J_RIGHT); - snprintf(tmp, 31, "%i", f_max); - show_text_simple(cr, tmp, -10.0, 0.0, NULL, 0.0, J_RIGHT); - - show_text_simple(cr, "0.00", 0.0, g_height+10.0, - NULL, -M_PI/3.0, J_RIGHT); - snprintf(tmp, 32, "%5.2f", osf_max); - show_text_simple(cr, tmp, g_width, g_height+10.0, - NULL, -M_PI/3.0, J_RIGHT); - - for ( b=0; b<nbins; b++ ) { - - double bar_height; - - bar_height = ((double)counts[b]/f_max)*g_height; - - cairo_new_path(cr); - cairo_rectangle(cr, (g_width/nbins)*b, g_height, - g_width/nbins, -bar_height); - cairo_set_source_rgb(cr, 0.0, 0.0, 1.0); - cairo_set_line_width(cr, 1.0); - cairo_stroke(cr); - - } - - cairo_new_path(cr); - cairo_rectangle(cr, 0.0, 0.0, g_width, g_height); - cairo_set_source_rgb(cr, 0.0, 0.0, 0.0); - cairo_set_line_width(cr, 1.0); - cairo_stroke(cr); -} - - -static void intensity_histogram(cairo_t *cr, Crystal **crystals, int n, - RefList *full, - signed int h, signed int k, signed int l) -{ - int f_max; - int i, b; - const int nbins = 30; - double int_max, int_inc; - double int_low[nbins]; - double int_high[nbins]; - int counts[nbins]; - const double g_width = 115.0; - const double g_height = 55.0; - char tmp[64]; - Reflection *f; - double Ifull, dsd_Ifull, pos, mI, bit; - int have_full; - - f = find_refl(full, h, k, l); - if ( f != NULL ) { - Ifull = get_intensity(f); - dsd_Ifull = get_esd_intensity(f) * sqrt(get_redundancy(f)); - have_full = 1; - } else { - Ifull = 0.0; - dsd_Ifull = 0.0; - have_full = 0; - } - - snprintf(tmp, 63, "%i %i %i", h, k, l); - show_text_simple(cr, tmp, g_width/2.0, -10.0, - "Sans Bold 10", 0.0, J_CENTER); - - int_max = 0.0; - int nmeas = 0; - for ( i=0; i<n; i++ ) { - - double osf; - Crystal *cryst = crystals[i]; - - if ( crystal_get_user_flag(cryst) ) continue; - - osf = crystal_get_osf(cryst); - - for ( f = find_refl(crystal_get_reflections(cryst), h, k, l); - f != NULL; - f = next_found_refl(f) ) - { - double Iobs, pcalc, Ifull_est; - - pcalc = get_partiality(f); - Iobs = get_intensity(f); - Ifull_est = Iobs / (pcalc * osf); - - if ( Ifull_est > int_max ) int_max = Ifull_est; - nmeas++; - } - - } - int_max *= 1.1; - int_inc = int_max / nbins; - - for ( b=0; b<nbins; b++ ) { - int_low[b] = b*int_inc; - int_high[b] = (b+1)*int_inc; - counts[b] = 0; - } - - for ( i=0; i<n; i++ ) { - - double osf; - Crystal *cryst = crystals[i]; - - if ( crystal_get_user_flag(cryst) ) continue; - - osf = crystal_get_osf(cryst); - - for ( f = find_refl(crystal_get_reflections(cryst), h, k, l); - f != NULL; - f = next_found_refl(f) ) - { - double Iobs, pcalc, Ifull_est; - - pcalc = get_partiality(f); - Iobs = get_intensity(f); - Ifull_est = Iobs / (pcalc * osf); - - for ( b=0; b<nbins; b++ ) { - if ( (Ifull_est >= int_low[b]) - && (Ifull_est < int_high[b]) ) { - counts[b]++; - break; - } - } - } - - } - - f_max = 0; - for ( b=0; b<nbins; b++ ) { - if ( counts[b] > f_max ) f_max = counts[b]; - } - f_max = (f_max/10)*10 + 10; - - snprintf(tmp, 32, "Max I=%.0f", int_max); - show_text_simple(cr, tmp, 10.0, 10.0, "Sans 9", 0.0, J_LEFT); - - for ( b=0; b<nbins; b++ ) { - - double bar_height; - - bar_height = ((double)counts[b]/f_max)*g_height; - - cairo_new_path(cr); - cairo_rectangle(cr, (g_width/nbins)*b, g_height, - g_width/nbins, -bar_height); - cairo_set_source_rgb(cr, 0.0, 0.0, 1.0); - cairo_set_line_width(cr, 1.0); - cairo_stroke(cr); - - } - - cairo_new_path(cr); - cairo_rectangle(cr, 0.0, 0.0, g_width, g_height); - cairo_set_source_rgb(cr, 0.0, 0.0, 0.0); - cairo_set_line_width(cr, 1.5); - cairo_stroke(cr); - - mI = int_high[nbins-1]; - pos = Ifull/mI; - bit = g_height / 15.0; - cairo_arc(cr, g_width*pos, g_height+bit, bit, 0.0, 2.0*M_PI); - if ( have_full ) { - cairo_set_source_rgb(cr, 0.0, 0.67, 0.45); - } else { - cairo_set_source_rgb(cr, 0.86, 0.0, 0.0); - } - cairo_fill(cr); - - if ( have_full ) { - - double eW = g_width*dsd_Ifull/mI; - - cairo_new_path(cr); - cairo_rectangle(cr, 0.0, g_height+bit*2.0, - g_width, g_height+bit*2.0); - //cairo_clip(cr); - - cairo_new_path(cr); - cairo_move_to(cr, g_width*pos - eW, g_height+bit); - cairo_line_to(cr, g_width*pos + eW, g_height+bit); - cairo_set_source_rgb(cr, 0.0, 0.67, 0.45); - cairo_set_line_width(cr, 2.0); - cairo_stroke(cr); - - cairo_reset_clip(cr); - } -} - - -static void watermark(struct _srcontext *sr) -{ - show_text(sr->cr, "Written by partialator from CrystFEL" - " version "PACKAGE_VERSION, sr->h-15.0, J_RIGHT, - "Sans 7"); -} - - -static void new_page(struct _srcontext *sr) -{ - cairo_show_page(sr->cr); - watermark(sr); -} - - -static void find_most_sampled_reflections(RefList *list, int n, signed int *h, - signed int *k, signed int *l) -{ - Reflection *refl; - RefListIterator *iter; - int *samples; - int j; - - for ( j=0; j<n; j++ ) { - h[j] = 0; - k[j] = 0; - l[j] = 0; - } - - samples = calloc(n, sizeof(int)); - - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - int red; - int i; - - red = get_redundancy(refl); - - for ( i=0; i<n; i++ ) { - - if ( red > samples[i] ) { - - int j; - - /* Shift everything down */ - for ( j=n-2; j>i; j-- ) { - h[j+1] = h[j]; - k[j+1] = k[j]; - l[j+1] = l[j]; - samples[j+1] = samples[j]; - } - - /* Add this in its place */ - get_indices(refl, &h[i], &k[i], &l[i]); - samples[i] = red; - - /* Don't compare against the others */ - break; - - } - - } - - } - - free(samples); -} - - - -SRContext *sr_titlepage(Crystal **crystals, int n, - const char *filename, const char *stream_filename, - const char *cmdline) -{ - char tmp[1024]; - struct _srcontext *sr; - - sr = malloc(sizeof(*sr)); - if ( sr == NULL ) return NULL; - - sr->w = PAGE_WIDTH; - sr->h = 595.0; - - sr->surf = cairo_pdf_surface_create(filename, sr->w, sr->h); - - if ( cairo_surface_status(sr->surf) != CAIRO_STATUS_SUCCESS ) { - fprintf(stderr, "Couldn't create Cairo surface\n"); - cairo_surface_destroy(sr->surf); - free(sr); - return NULL; - } - - sr->cr = cairo_create(sr->surf); - watermark(sr); - - snprintf(tmp, 1023, "%s", stream_filename); - show_text(sr->cr, tmp, 10.0, J_CENTER, "Sans Bold 16"); - snprintf(tmp, 1023, "partialator %s", cmdline); - show_text(sr->cr, tmp, 45.0, J_LEFT, "Mono 7"); - - return sr; -} - - -void sr_iteration(SRContext *sr, int iteration, struct srdata *d) -{ - int i; - char page_title[1024]; - double dash[] = {2.0, 2.0}; - - if ( sr == NULL ) return; - - snprintf(page_title, 1023, "After %i iteration%s", - iteration, iteration==1?"":"s"); - - new_page(sr); - show_text(sr->cr, page_title, 10.0, J_CENTER, "Sans Bold 16"); - - cairo_save(sr->cr); - cairo_translate(sr->cr, 480.0, 350.0); - scale_factor_histogram(sr->cr, d->crystals, d->n, - "Distribution of overall scale factors"); - cairo_restore(sr->cr); - - /* Draw partiality plots (three graphs together) */ - cairo_save(sr->cr); - - cairo_translate(sr->cr, 70.0, 330.0); - partiality_graph(sr->cr, d->crystals, d->n, d->full); - - cairo_save(sr->cr); - cairo_move_to(sr->cr, 0.0, 0.0); - cairo_line_to(sr->cr, 0.0, -30.0); - cairo_move_to(sr->cr, 200.0, 0.0); - cairo_line_to(sr->cr, 200.0, -30.0); - cairo_set_dash(sr->cr, dash, 2, 0.0); - cairo_stroke(sr->cr); - cairo_set_dash(sr->cr, NULL, 0, 0.0); - cairo_translate(sr->cr, 0.0, -150.0); - partiality_histogram(sr->cr, d->crystals, d->n, d->full, 1, 0); - cairo_restore(sr->cr); - - cairo_save(sr->cr); - cairo_move_to(sr->cr, 200.0, 0.0); - cairo_line_to(sr->cr, 230.0, 00.0); - cairo_move_to(sr->cr, 200.0, 200.0); - cairo_line_to(sr->cr, 230.0, 200.0); - cairo_set_dash(sr->cr, dash, 2, 0.0); - cairo_stroke(sr->cr); - cairo_set_dash(sr->cr, NULL, 0, 0.0); - cairo_translate(sr->cr, 230.0, 200.0); - cairo_rotate(sr->cr, -M_PI_2); - partiality_histogram(sr->cr, d->crystals, d->n, d->full, 0, 1); - cairo_restore(sr->cr); - - cairo_restore(sr->cr); - - if ( iteration == 0 ) { - find_most_sampled_reflections(d->full, 9, - sr->ms_h, sr->ms_k, sr->ms_l); - } - - for ( i=0; i<9; i++ ) { - - int x, y; - - x = i % 3; - y = i / 3; - - cairo_save(sr->cr); - cairo_translate(sr->cr, 400.0+140.0*x, 60.0+80.0*y); - intensity_histogram(sr->cr, d->crystals, d->n, d->full, - sr->ms_h[i], sr->ms_k[i], sr->ms_l[i]); - cairo_restore(sr->cr); - - } -} - - -void sr_finish(SRContext *sr) -{ - cairo_surface_finish(sr->surf); - cairo_destroy(sr->cr); -} diff --git a/src/scaling-report.h b/src/scaling-report.h deleted file mode 100644 index 55037038..00000000 --- a/src/scaling-report.h +++ /dev/null @@ -1,82 +0,0 @@ -/* - * scaling-report.h - * - * Write a nice PDF of scaling parameters - * - * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2011-2014 Thomas White <taw@physics.org> - * - * 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 SCALING_REPORT_H -#define SCALING_REPORT_H - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - - -#include "utils.h" - -typedef struct _srcontext SRContext; /* Opaque */ - -/* Information is logged in this structure */ -struct srdata -{ - Crystal **crystals; - int n; - RefList *full; - - int n_filtered; - int n_refined; -}; - -#if defined(HAVE_CAIRO) && defined(HAVE_PANGO) && defined(HAVE_PANGOCAIRO) - -extern SRContext *sr_titlepage(Crystal **crystals, int n, - const char *filename, - const char *stream_filename, - const char *cmdline); - -extern void sr_iteration(SRContext *sr, int iteration, struct srdata *d); - -extern void sr_finish(SRContext *sr); - -#else /* defined(HAVE_CAIRO) && defined(HAVE_PANGO) && ... */ - -SRContext *sr_titlepage(Crystal **crystals, int n, const char *filename, - const char *stream_filename, const char *cmdline) -{ - return NULL; -} - -void sr_iteration(SRContext *sr, int iteration, struct srdata *d) -{ -} - -void sr_finish(SRContext *sr) -{ -} - -#endif /* defined(HAVE_CAIRO) && defined(HAVE_PANGO) && ... */ - - -#endif /* SCALING_REPORT_H */ |