diff options
Diffstat (limited to 'libcrystfel/src/predict-refine.c')
-rw-r--r-- | libcrystfel/src/predict-refine.c | 484 |
1 files changed, 361 insertions, 123 deletions
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 19c689cd..43b54092 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -3,11 +3,11 @@ * * Prediction refinement * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2023 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2020 Thomas White <taw@physics.org> + * 2010-2023 Thomas White <taw@physics.org> * 2016 Valerio Mariani * * This file is part of CrystFEL. @@ -33,91 +33,325 @@ #include <assert.h> #include <gsl/gsl_matrix.h> #include <gsl/gsl_vector.h> +#include <gsl/gsl_linalg.h> #include "image.h" #include "geometry.h" #include "cell-utils.h" +#include "predict-refine.h" +#include "profile.h" +#include "crystfel-mille.h" /** \file predict-refine.h */ +double r_dev(struct reflpeak *rp) +{ + /* Excitation error term */ + return get_exerr(rp->refl); +} -/* 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) +double fs_dev(struct reflpeak *rp, struct detgeom *det) +{ + double fsh, ssh; + get_detector_pos(rp->refl, &fsh, &ssh); + return fsh - rp->peak->fs; +} + -/* Parameters to refine */ -static const enum gparam rv[] = +double ss_dev(struct reflpeak *rp, struct detgeom *det) { - 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 detgeom_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 detgeom_panel *p, double dx, double dy) + double fsh, ssh; + get_detector_pos(rp->refl, &fsh, &ssh); + return ssh - rp->peak->ss; +} + + +double r_gradient(int param, Reflection *refl, UnitCell *cell, double wavelength) { - double xs, ys; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double xl, yl, zl; + signed int hs, ks, ls; + double tl, phi, azi; + + get_symmetric_indices(refl, &hs, &ks, &ls); + + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + xl = hs*asx + ks*bsx + ls*csx; + yl = hs*asy + ks*bsy + ls*csy; + zl = hs*asz + ks*bsz + ls*csz; + + tl = sqrt(xl*xl + yl*yl); + phi = angle_between_2d(tl, zl+1.0/wavelength, 0.0, 1.0); /* 2theta */ + azi = atan2(yl, xl); /* azimuth */ + + switch ( param ) { + + case GPARAM_ASX : + return - hs * sin(phi) * cos(azi); + + case GPARAM_BSX : + return - ks * sin(phi) * cos(azi); + + case GPARAM_CSX : + return - ls * sin(phi) * cos(azi); + + case GPARAM_ASY : + return - hs * sin(phi) * sin(azi); - xs = fs*p->fsx + ss*p->ssx; /* pixels */ - ys = fs*p->fsy + ss*p->ssy; /* pixels */ + case GPARAM_BSY : + return - ks * sin(phi) * sin(azi); - *px = (xs + p->cnx) * p->pixel_pitch + dx; /* metres */ - *py = (ys + p->cny) * p->pixel_pitch + dy; /* metres */ + case GPARAM_CSY : + return - ls * sin(phi) * sin(azi); + + case GPARAM_ASZ : + return - hs * cos(phi); + + case GPARAM_BSZ : + return - ks * cos(phi); + + case GPARAM_CSZ : + return - ls * cos(phi); + + /* Detector movements don't affect excitation error */ + case GPARAM_DET_TX : + case GPARAM_DET_TY : + case GPARAM_DET_TZ : + case GPARAM_DET_RX : + case GPARAM_DET_RY : + case GPARAM_DET_RZ : + return 0.0; + + } + + ERROR("No r gradient defined for parameter %i\n", param); + abort(); } -static double r_dev(struct reflpeak *rp) +/* Spot position gradients for diffraction physics (anything that changes the + * diffracted ray direction) */ +int fs_ss_gradient_physics(int param, Reflection *refl, UnitCell *cell, + struct detgeom_panel *p, gsl_matrix *Minv, + double fs, double ss, double mu, + float *fsg, float *ssg) { - /* Excitation error term */ - return get_exerr(rp->refl); + signed int h, k, l; + gsl_vector *dRdp; + gsl_vector *v; + + get_indices(refl, &h, &k, &l); + + dRdp = gsl_vector_calloc(3); + + switch ( param ) { + + case GPARAM_ASX : + gsl_vector_set(dRdp, 0, h); + break; + + case GPARAM_BSX : + gsl_vector_set(dRdp, 0, k); + break; + + case GPARAM_CSX : + gsl_vector_set(dRdp, 0, l); + break; + + case GPARAM_ASY : + gsl_vector_set(dRdp, 1, h); + break; + + case GPARAM_BSY : + gsl_vector_set(dRdp, 1, k); + break; + + case GPARAM_CSY : + gsl_vector_set(dRdp, 1, l); + break; + + case GPARAM_ASZ : + gsl_vector_set(dRdp, 2, h); + break; + + case GPARAM_BSZ : + gsl_vector_set(dRdp, 2, k); + break; + + case GPARAM_CSZ : + gsl_vector_set(dRdp, 2, l); + break; + + default : + ERROR("Invalid physics gradient %i\n", param); + return 1; + } + + v = gsl_vector_calloc(3); + gsl_blas_dgemv(CblasNoTrans, 1.0, Minv, dRdp, 0.0, v); + + *fsg = mu*(gsl_vector_get(v, 1) - fs*gsl_vector_get(v, 0)); + *ssg = mu*(gsl_vector_get(v, 2) - ss*gsl_vector_get(v, 0)); + + gsl_vector_free(v); + gsl_vector_free(dRdp); + + return 0; } -static double x_dev(struct reflpeak *rp, struct detgeom *det, - double dx, double dy) +/* Spot position gradients for panel motions (translation or rotation) */ +int fs_ss_gradient_panel(int param, Reflection *refl, UnitCell *cell, + struct detgeom_panel *p, gsl_matrix *Minv, + double fs, double ss, double mu, + gsl_vector *t, double cx, double cy, double cz, + float *fsg, float *ssg) { - /* Peak position term */ - double xpk, ypk, xh, yh; - double fsh, ssh; - twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel, dx, dy); - get_detector_pos(rp->refl, &fsh, &ssh); - twod_mapping(fsh, ssh, &xh, &yh, rp->panel, dx, dy); - return xh-xpk; + gsl_vector *v; + gsl_matrix *gM; /* M^-1 * dM/dx * M^-1 */ + gsl_matrix *dMdp = gsl_matrix_calloc(3, 3); + + switch ( param ) { + + case GPARAM_DET_TX : + gsl_matrix_set(dMdp, 0, 0, 1.0); + break; + + case GPARAM_DET_TY : + gsl_matrix_set(dMdp, 1, 0, 1.0); + break; + + case GPARAM_DET_TZ : + gsl_matrix_set(dMdp, 2, 0, 1.0); + break; + + case GPARAM_DET_RX : + gsl_matrix_set(dMdp, 1, 0, cz-p->pixel_pitch*p->cnz); + gsl_matrix_set(dMdp, 2, 0, p->pixel_pitch*p->cny-cy); + gsl_matrix_set(dMdp, 1, 1, -p->pixel_pitch*p->fsz); + gsl_matrix_set(dMdp, 2, 1, p->pixel_pitch*p->fsy); + gsl_matrix_set(dMdp, 1, 2, -p->pixel_pitch*p->ssz); + gsl_matrix_set(dMdp, 2, 2, p->pixel_pitch*p->ssy); + break; + + case GPARAM_DET_RY : + gsl_matrix_set(dMdp, 0, 0, p->pixel_pitch*p->cnz-cz); + gsl_matrix_set(dMdp, 2, 0, cx-p->pixel_pitch*p->cnx); + gsl_matrix_set(dMdp, 0, 1, p->pixel_pitch*p->fsz); + gsl_matrix_set(dMdp, 2, 1, -p->pixel_pitch*p->fsx); + gsl_matrix_set(dMdp, 0, 2, p->pixel_pitch*p->ssz); + gsl_matrix_set(dMdp, 2, 2, -p->pixel_pitch*p->ssx); + break; + + case GPARAM_DET_RZ : + gsl_matrix_set(dMdp, 0, 0, cy-p->pixel_pitch*p->cny); + gsl_matrix_set(dMdp, 1, 0, p->pixel_pitch*p->cnx-cx); + gsl_matrix_set(dMdp, 0, 1, -p->pixel_pitch*p->fsy); + gsl_matrix_set(dMdp, 1, 1, p->pixel_pitch*p->fsx); + gsl_matrix_set(dMdp, 0, 2, -p->pixel_pitch*p->ssy); + gsl_matrix_set(dMdp, 1, 2, p->pixel_pitch*p->ssx); + break; + + default: + ERROR("Invalid panel gradient %i\n", param); + return 1; + + } + + gM = matrix_mult3(Minv, dMdp, Minv); + gsl_matrix_free(dMdp); + + v = gsl_vector_calloc(3); + gsl_blas_dgemv(CblasNoTrans, -1.0, gM, t, 0.0, v); + gsl_vector_free(t); + gsl_matrix_free(gM); + + *fsg = mu*(gsl_vector_get(v, 1) - fs*gsl_vector_get(v, 0)); + *ssg = mu*(gsl_vector_get(v, 2) - ss*gsl_vector_get(v, 0)); + + gsl_vector_free(v); + + return 0; } -static double y_dev(struct reflpeak *rp, struct detgeom *det, - double dx, double dy) +/* Returns the gradient of fs_dev and ss_dev w.r.t. any parameter. + * cx,cy,cz are the rotation axis coordinates (only 2 in use at any time) + * in metres (not pixels) */ +int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell, + struct detgeom_panel *p, gsl_matrix *Minv, + double cx, double cy, double cz, + float *fsg, float *ssg) { - /* Peak position term */ - double xpk, ypk, xh, yh; - double fsh, ssh; - twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel, dx, dy); - get_detector_pos(rp->refl, &fsh, &ssh); - twod_mapping(fsh, ssh, &xh, &yh, rp->panel, dx, dy); - return yh-ypk; + signed int h, k, l; + double xl, yl, zl, kpred; + double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; + gsl_vector *t; + gsl_vector *v; + gsl_matrix *M; + double mu; + double fs, ss; + + get_indices(refl, &h, &k, &l); + kpred = get_kpred(refl); + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + xl = h*asx + k*bsx + l*csx; + yl = h*asy + k*bsy + l*csy; + zl = h*asz + k*bsz + l*csz; + + /* Set up matrix equation */ + M = gsl_matrix_alloc(3, 3); + v = gsl_vector_alloc(3); + t = gsl_vector_alloc(3); + if ( (M==NULL) || (v==NULL) || (t==NULL) ) { + ERROR("Failed to allocate vectors for gradient calculation\n"); + return 1; + } + + gsl_vector_set(t, 0, xl); + gsl_vector_set(t, 1, yl); + gsl_vector_set(t, 2, kpred+zl); + + gsl_matrix_set(M, 0, 0, (p->cnx)*p->pixel_pitch); + gsl_matrix_set(M, 0, 1, (p->fsx)*p->pixel_pitch); + gsl_matrix_set(M, 0, 2, (p->ssx)*p->pixel_pitch); + gsl_matrix_set(M, 1, 0, (p->cny)*p->pixel_pitch); + gsl_matrix_set(M, 1, 1, (p->fsy)*p->pixel_pitch); + gsl_matrix_set(M, 1, 2, (p->ssy)*p->pixel_pitch); + gsl_matrix_set(M, 2, 0, (p->cnz)*p->pixel_pitch); + gsl_matrix_set(M, 2, 1, (p->fsz)*p->pixel_pitch); + gsl_matrix_set(M, 2, 2, (p->ssz)*p->pixel_pitch); + + if ( gsl_linalg_HH_solve(M, t, v) ) { + ERROR("Failed to solve gradient equation\n"); + return 1; + } + gsl_matrix_free(M); + + mu = 1.0 / gsl_vector_get(v, 0); + fs = mu*gsl_vector_get(v, 1); + ss = mu*gsl_vector_get(v, 2); + gsl_vector_free(v); + + if ( param <= GPARAM_CSZ ) { + gsl_vector_free(t); + return fs_ss_gradient_physics(param, refl, cell, p, + Minv, fs, ss, mu, + fsg, ssg); + } else { + return fs_ss_gradient_panel(param, refl, cell, p, + Minv, fs, ss, mu, t, + cx, cy, cz, + fsg, ssg); + } } @@ -244,7 +478,6 @@ static int pair_peaks(struct image *image, Crystal *cr, rps[n].refl = refl; rps[n].peak = f; - rps[n].panel = &image->detgeom->panels[f->pn]; n++; } @@ -349,8 +582,7 @@ int refine_radius(Crystal *cr, struct image *image) static int iterate(struct reflpeak *rps, int n, UnitCell *cell, - struct image *image, - double *total_x, double *total_y, double *total_z) + struct image *image, gsl_matrix **Minvs) { int i; gsl_matrix *M; @@ -359,6 +591,18 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; + const enum gparam rv[] = { + GPARAM_ASX, + GPARAM_ASY, + GPARAM_ASZ, + GPARAM_BSX, + GPARAM_BSY, + GPARAM_BSZ, + GPARAM_CSX, + GPARAM_CSY, + GPARAM_CSZ, + }; + const int num_params = 9; /* Number of parameters to refine */ M = gsl_matrix_calloc(num_params, num_params); @@ -367,17 +611,21 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, for ( i=0; i<n; i++ ) { int k; - double gradients[num_params]; + float fs_gradients[num_params]; + float ss_gradients[num_params]; + float r_gradients[num_params]; double w; - /* Excitation error terms */ - w = EXC_WEIGHT * rps[i].Ih; - + /* Calculate all gradients for this parameter */ for ( k=0; k<num_params; k++ ) { - gradients[k] = r_gradient(cell, rv[k], rps[i].refl, - image); + int pn = rps[i].peak->pn; + r_gradients[k] = r_gradient(rv[k], rps[i].refl, cell, image->lambda); + fs_ss_gradient(rv[k], rps[i].refl, cell, &image->detgeom->panels[pn], + Minvs[pn], 0, 0, 0, &fs_gradients[k], &ss_gradients[k]); } + /* Excitation error terms */ + w = EXC_WEIGHT * rps[i].Ih; for ( k=0; k<num_params; k++ ) { int g; @@ -390,7 +638,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, /* Matrix is symmetric */ if ( g > k ) continue; - M_c = w * gradients[g] * gradients[k]; + M_c = w * r_gradients[g] * r_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); @@ -398,18 +646,13 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, } v_c = w * r_dev(&rps[i]); - v_c *= -gradients[k]; + v_c *= -r_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].refl, cell, - rps[i].panel); - } - + /* Positional fs terms */ for ( k=0; k<num_params; k++ ) { int g; @@ -422,26 +665,21 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, /* Matrix is symmetric */ if ( g > k ) continue; - M_c = gradients[g] * gradients[k]; + M_c = fs_gradients[g] * fs_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->detgeom, *total_x, *total_y); - v_c *= -gradients[k]; + v_c = fs_dev(&rps[i], image->detgeom); + v_c *= -fs_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].refl, cell, - rps[i].panel); - } - + /* Positional ss terms */ for ( k=0; k<num_params; k++ ) { int g; @@ -454,15 +692,15 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, /* Matrix is symmetric */ if ( g > k ) continue; - M_c = gradients[g] * gradients[k]; + M_c = ss_gradients[g] * ss_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->detgeom, *total_x, *total_y); - v_c *= -gradients[k]; + v_c = ss_dev(&rps[i], image->detgeom); + v_c *= -ss_gradients[k]; v_curr = gsl_vector_get(v, k); gsl_vector_set(v, k, v_curr + v_c); @@ -472,14 +710,8 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, int k; for ( k=0; k<num_params; k++ ) { - double M_curr; - M_curr = gsl_matrix_get(M, k, k); - if ( (rv[k] == GPARAM_DETX) || (rv[k] == GPARAM_DETY) ) { - M_curr += 10.0; - } else { - M_curr += 1e-18; - } - gsl_matrix_set(M, k, k, M_curr); + double M_curr = gsl_matrix_get(M, k, k); + gsl_matrix_set(M, k, k, M_curr+1e-7); } //show_matrix_eqn(M, v); @@ -513,9 +745,6 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, csx += gsl_vector_get(shifts, 6); csy += gsl_vector_get(shifts, 7); csz += gsl_vector_get(shifts, 8); - *total_x += gsl_vector_get(shifts, 9); - *total_y += gsl_vector_get(shifts, 10); - *total_z += 0.0; cell_set_reciprocal(cell, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); @@ -527,8 +756,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell, } -static double pred_residual(struct reflpeak *rps, int n, struct detgeom *det, - double dx, double dy) +static double pred_residual(struct reflpeak *rps, int n, struct detgeom *det) { int i; double res = 0.0; @@ -542,13 +770,13 @@ static double pred_residual(struct reflpeak *rps, int n, struct detgeom *det, r = 0.0; for ( i=0; i<n; i++ ) { - r += pow(x_dev(&rps[i], det, dx, dy), 2.0); + r += pow(det->panels[rps[i].peak->pn].pixel_pitch*fs_dev(&rps[i], det), 2.0); } res += r; r = 0.0; for ( i=0; i<n; i++ ) { - r += pow(y_dev(&rps[i], det, dx, dy), 2.0); + r += pow(det->panels[rps[i].peak->pn].pixel_pitch*ss_dev(&rps[i], det), 2.0); } res += r; @@ -569,18 +797,15 @@ static void free_rps_noreflist(struct reflpeak *rps, int n) } -int refine_prediction(struct image *image, Crystal *cr) +int refine_prediction(struct image *image, Crystal *cr, Mille *mille) { 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; - double orig_shift_x, orig_shift_y; char tmp[256]; + gsl_matrix **Minvs; rps = malloc(image_feature_count(image->features) * sizeof(struct reflpeak)); @@ -595,9 +820,7 @@ int refine_prediction(struct image *image, Crystal *cr) } crystal_set_reflections(cr, reflist); - crystal_get_det_shift(cr, &total_x, &total_y); - orig_shift_x = total_x; - orig_shift_y = total_y; + Minvs = make_panel_minvs(image->detgeom); /* Normalise the intensities to max 1 */ max_I = -INFINITY; @@ -620,29 +843,36 @@ int refine_prediction(struct image *image, Crystal *cr) } //STATUS("Initial residual = %e\n", - // pred_residual(rps, n, image->detgeom, total_x, total_y)); + // pred_residual(rps, n, image->detgeom)); - /* Refine */ - for ( i=0; i<MAX_CYCLES; i++ ) { + /* Refine (max 10 cycles) */ + for ( i=0; i<10; i++ ) { update_predictions(cr); - if ( iterate(rps, n, crystal_get_cell(cr), image, - &total_x, &total_y, &total_z) ) + if ( iterate(rps, n, crystal_get_cell(cr), image, Minvs) ) { crystal_set_reflections(cr, NULL); return 1; } - crystal_set_det_shift(cr, total_x, total_y); //STATUS("Residual after %i = %e\n", i, - // pred_residual(rps, n, image->detgeom, total_x, total_y)); + // pred_residual(rps, n, image->detgeom)); } //STATUS("Final residual = %e\n", - // pred_residual(rps, n, image->detgeom, total_x, total_y)); + // pred_residual(rps, n, image->detgeom)); snprintf(tmp, 255, "predict_refine/final_residual = %e", - pred_residual(rps, n, image->detgeom, total_x, total_y)); + pred_residual(rps, n, image->detgeom)); crystal_add_notes(cr, tmp); - crystal_set_det_shift(cr, total_x, total_y); + if ( mille != NULL ) { + profile_start("mille-calc"); + write_mille(mille, n, crystal_get_cell(cr), rps, image, Minvs); + profile_end("mille-calc"); + } + + for ( i=0; i<image->detgeom->n_panels; i++ ) { + gsl_matrix_free(Minvs[i]); + } + free(Minvs); crystal_set_reflections(cr, NULL); reflist_free(reflist); @@ -650,9 +880,17 @@ int refine_prediction(struct image *image, Crystal *cr) n = pair_peaks(image, cr, NULL, rps); free_rps_noreflist(rps, n); if ( n < 10 ) { - crystal_set_det_shift(cr, orig_shift_x, orig_shift_y); + if ( mille != NULL ) { + crystfel_mille_delete_last_record(mille); + } return 1; } + if ( mille != NULL ) { + profile_start("mille-write"); + crystfel_mille_write_record(mille); + profile_end("mille-write"); + } + return 0; } |