From 0a7dcfd71ae2c267ac748c43637e27aef403db7e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 25 Jun 2015 17:56:45 +0200 Subject: Move predict-refine to libcrystfel --- libcrystfel/Makefile.am | 4 +- libcrystfel/src/predict-refine.c | 757 +++++++++++++++++++++++++++++++++++++++ libcrystfel/src/predict-refine.h | 45 +++ 3 files changed, 804 insertions(+), 2 deletions(-) create mode 100644 libcrystfel/src/predict-refine.c create mode 100644 libcrystfel/src/predict-refine.h (limited to 'libcrystfel') diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index 8e5acaaf..1a0c3efe 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -9,7 +9,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/reflist-utils.c src/filters.c \ src/render.c src/index.c src/dirax.c src/mosflm.c \ src/cell-utils.c src/integer_matrix.c src/crystal.c \ - src/xds.c src/integration.c \ + src/xds.c src/integration.c src/predict-refine.c \ src/histogram.c src/events.c if HAVE_FFTW @@ -28,7 +28,7 @@ libcrystfel_la_include_HEADERS = ${top_srcdir}/version.h \ src/filters.h src/dirax.h src/mosflm.h \ src/cell-utils.h \ src/integer_matrix.h src/crystal.h \ - src/xds.h \ + src/xds.h src/predict-refine.h \ src/integration.h src/histogram.h \ src/events.h src/asdf.h diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c new file mode 100644 index 00000000..2d845feb --- /dev/null +++ b/libcrystfel/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 + * + * 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 . + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include +#include + +#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; ifs; + 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; ifeatures); 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; ifs, 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; in_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 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; kdet, + image->lambda, cell); + } + + for ( k=0; k 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; kdet, + image->lambda, cell); + } + + for ( k=0; k 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; idet, 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; ifeatures) + * 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; iintensity; + 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; iintensity / max_I; + } + + //STATUS("Initial residual = %e\n", residual(rps, n, image->det)); + + /* Refine */ + for ( i=0; idet)); + } + //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/libcrystfel/src/predict-refine.h b/libcrystfel/src/predict-refine.h new file mode 100644 index 00000000..c763d2ca --- /dev/null +++ b/libcrystfel/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 + * + * 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 . + * + */ + +#ifndef PREDICT_REFINE_H +#define PREDICT_REFINE_H + + +#ifdef HAVE_CONFIG_H +#include +#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 */ -- cgit v1.2.3