diff options
author | Thomas White <taw@physics.org> | 2015-03-20 17:13:15 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2015-04-20 15:50:39 +0200 |
commit | e05288c1074be5928cca8581ab4d07d80853c22a (patch) | |
tree | 8241e81a13b6c06703ec53981bc56d0d29057c0e /src/predict-refine.c | |
parent | a97aeb4506de83a9e556caa54557179433c7687b (diff) |
refine_{radius,prediction}: unite common code (pair_peaks), avoid preliminary integration round
Diffstat (limited to 'src/predict-refine.c')
-rw-r--r-- | src/predict-refine.c | 166 |
1 files changed, 113 insertions, 53 deletions
diff --git a/src/predict-refine.c b/src/predict-refine.c index a2a4ce53..52496198 100644 --- a/src/predict-refine.c +++ b/src/predict-refine.c @@ -55,28 +55,77 @@ struct reflpeak { * (we assume this never changes) */ }; -static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist, - struct reflpeak *rps, struct detector *det) + +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; +} + + +/* 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; const double min_dist = 0.05; int n_acc = 0; - int n_notintegrated = 0; - for ( i=0; i<image_feature_count(flist); i++ ) { + for ( i=0; i<image_feature_count(image->features); i++ ) { struct imagefeature *f; double h, k, l, hd, kd, ld; /* Assume all image "features" are genuine peaks */ - f = image_get_feature(flist, i); + f = image_get_feature(image->features, i); if ( f == NULL ) continue; double ax, ay, az; double bx, by, bz; double cx, cy, cz; - cell_get_cartesian(cell, + cell_get_cartesian(crystal_get_cell(cr), &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); /* Decimal and fractional Miller indices of nearest @@ -94,17 +143,28 @@ static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist, && (fabs(l - ld) < min_dist) ) { Reflection *refl; + struct panel *p; - /* Dig out the reflection */ - refl = find_refl(reflist, h, k, l); + refl = reflection_new(h, k, l); if ( refl == NULL ) { - n_notintegrated++; - continue; + ERROR("Failed to create reflection\n"); + return 0; } + if ( reflist != NULL ) add_refl_to_list(refl, 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_acc].refl = refl; rps[n_acc].peak = f; - rps[n_acc].panel = find_panel(det, f->fs, f->ss); + rps[n_acc].panel = p; n_acc++; } @@ -114,19 +174,48 @@ static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist, } -static void twod_mapping(double fs, double ss, double *px, double *py, - struct panel *p) +static int cmpd2(const void *av, const void *bv) { - double xs, ys; + struct reflpeak *a, *b; - xs = fs*p->fsx + ss*p->ssx; - ys = fs*p->fsy + ss*p->ssy; + a = (struct reflpeak *)av; + b = (struct reflpeak *)bv; - *px = (xs + p->cnx) / p->res; - *py = (ys + p->cny) / p->res; + if ( fabs(r_dev(a)) < fabs(r_dev(b)) ) return -1; + return 1; } +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 + d(yh-ypk)/dP, where P = any parameter */ @@ -239,39 +328,6 @@ static double y_gradient(int param, struct reflpeak *rp, struct detector *det, } -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 int iterate(struct reflpeak *rps, int n, UnitCell *cell, struct image *image) { @@ -459,18 +515,20 @@ int refine_prediction(struct image *image, Crystal *cr) int i; struct reflpeak *rps; double max_I; + RefList *reflist; rps = malloc(image_feature_count(image->features) * sizeof(struct reflpeak)); if ( rps == NULL ) return 1; - n = pair_peaks(image->features, crystal_get_cell(cr), - crystal_get_reflections(cr), rps, image->det); + reflist = reflist_new(); + n = pair_peaks(image, cr, reflist, rps); if ( n < 10 ) { ERROR("Too few paired peaks to refine orientation.\n"); free(rps); return 1; } + crystal_set_reflections(cr, reflist); /* Normalise the intensities to max 1 */ max_I = -INFINITY; @@ -493,6 +551,8 @@ int refine_prediction(struct image *image, Crystal *cr) update_partialities(cr, PMODEL_SCSPHERE); } + crystal_set_reflections(cr, NULL); + reflist_free(reflist); free(rps); return 0; } |