aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-09-28 18:20:56 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:00 +0100
commit0dcf9c99fa6e3e78524b34ce24a40011c2554b3c (patch)
treeba4c0dea8d0e002590550b460fa1e63a74d65d7e
parent8f43ca33d5edd9f75660174870e16afe8b9913bd (diff)
Update sanity check
-rw-r--r--src/Makefile.am3
-rw-r--r--src/index.c11
-rw-r--r--src/index.h2
-rw-r--r--src/peaks.c63
-rw-r--r--src/reintegrate.c5
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");