diff options
| -rw-r--r-- | src/Makefile.am | 3 | ||||
| -rw-r--r-- | src/index.c | 11 | ||||
| -rw-r--r-- | src/index.h | 2 | ||||
| -rw-r--r-- | src/peaks.c | 63 | ||||
| -rw-r--r-- | src/reintegrate.c | 5 |
5 files changed, 68 insertions, 16 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index aba2d6fc..e147b78e 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -67,7 +67,8 @@ cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c diffraction.c \ cubeit_LDADD = @LIBS@ reintegrate_SOURCES = reintegrate.c cell.c hdf5-file.c utils.c detector.c \ - peaks.c image.c diffraction.c sfac.c geometry.c stream.c + peaks.c image.c diffraction.c sfac.c geometry.c stream.c \ + index.c dirax.c templates.c symmetry.c reintegrate_LDADD = @LIBS@ INCLUDES = "-I$(top_srcdir)/data" diff --git a/src/index.c b/src/index.c index b5d9da9e..8d4b5dd6 100644 --- a/src/index.c +++ b/src/index.c @@ -99,8 +99,7 @@ static void write_drx(struct image *image) } -void index_pattern(struct image *image, UnitCell *cell, IndexingMethod indm, - int no_match, int verbose, IndexingPrivate *ipriv) +void map_all_peaks(struct image *image) { int i; int nc = 0; @@ -121,7 +120,15 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod indm, if ( nc ) { ERROR("Failed to map %i reflections\n", nc); } +} + + +void index_pattern(struct image *image, UnitCell *cell, IndexingMethod indm, + int no_match, int verbose, IndexingPrivate *ipriv) +{ + int i; + map_all_peaks(image); write_drx(image); image->ncells = 0; diff --git a/src/index.h b/src/index.h index 2e957538..259d3eb6 100644 --- a/src/index.h +++ b/src/index.h @@ -35,6 +35,8 @@ extern IndexingPrivate *prepare_indexing(IndexingMethod indm, UnitCell *cell, const char *filename, struct detector *det); +extern void map_all_peaks(struct image *image); + extern void index_pattern(struct image *image, UnitCell *cell, IndexingMethod indm, int no_match, int verbose, IndexingPrivate *priv); diff --git a/src/peaks.c b/src/peaks.c index 16a5f2a2..cd414782 100644 --- a/src/peaks.c +++ b/src/peaks.c @@ -549,30 +549,69 @@ int peak_sanity_check(struct image *image, UnitCell *cell, int circular_domain, double domain_r) { int i; + int n_feat = 0; int n_sane = 0; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + double alen, blen, clen; find_projected_peaks(image, cell, circular_domain, domain_r); if ( image->n_hits == 0 ) return 0; /* Failed sanity check: no peaks */ - for ( i=0; i<image->n_hits; i++ ) { + /* "Borrow" direction values to get reciprocal lengths */ + cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + alen = modulus(ax, ay, az); + blen = modulus(bx, by, bz); + clen = modulus(cx, cy, cz); - double d; - int idx; + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + fesetround(1); /* Round towards nearest */ + n_feat = image_feature_count(image->features); + for ( i=0; i<n_feat; i++ ) { + + double dist; + struct rvec q; struct imagefeature *f; + double hd, kd, ld; + signed int h, k, l; + double dh, dk, dl; + + f = image_get_feature(image->features, i); + + /* Get closest hkl */ + q = get_q(image, f->x, f->y, 1, NULL, 1.0/image->lambda); + + hd = q.u * ax + q.v * ay + q.w * az; + kd = q.u * bx + q.v * by + q.w * bz; + ld = q.u * cx + q.v * cy + q.w * cz; + + h = lrint(hd); k = lrint(kd); l = lrint(ld); + + dh = hd - h; dk = kd - k; dl = ld - l; + + if ( circular_domain ) { - f = image_feature_closest(image->features, - image->hits[i].x, image->hits[i].y, - &d, &idx); - if ( (f != NULL) && (d < PEAK_CLOSE) ) { - n_sane++; + /* Circular integration domain */ + dist = sqrt(pow(dh*alen, 2.0) + pow(dk*blen, 2.0) + + pow(dl*clen, 2.0)); + if ( dist <= domain_r ) n_sane++; + + } else { + + /* "Crystallographic" integration domain */ + dist = sqrt(pow(dh, 2.0) + pow(dk, 2.0) + pow(dl, 2.0)); + if ( dist <= domain_r ) n_sane++; } + + } - STATUS("Sanity factor: %f / %f = %f\n", (float)n_sane, - (float)image->n_hits, - (float)n_sane / (float)image->n_hits); - if ( (float)n_sane / (float)image->n_hits < 0.8 ) return 0; + STATUS("Sanity factor: %f / %f = %f\n", (float)n_sane, (float)n_feat, + (float)n_sane / (float)n_feat); + if ( (float)n_sane / (float)n_feat < 0.8 ) return 0; return 1; } diff --git a/src/reintegrate.c b/src/reintegrate.c index 53117163..a278a3ec 100644 --- a/src/reintegrate.c +++ b/src/reintegrate.c @@ -29,6 +29,7 @@ #include "symmetry.h" #include "peaks.h" #include "stream.h" +#include "index.h" #define MAX_THREADS (256) @@ -121,9 +122,11 @@ static void process_image(struct process_args *pargs) hdf5_read(hdfile, &image, pargs->config_satcorr); + map_all_peaks(&image); + /* Sanity check */ if ( pargs->config_sanity - && !peak_sanity_check(&image, image.indexed_cell, 0, 0.1) ) { + && !peak_sanity_check(&image, image.indexed_cell, 1, 0.006e9) ) { STATUS("Failed peak sanity check.\n"); |