diff options
| -rw-r--r-- | src/templates.c | 81 |
1 files changed, 78 insertions, 3 deletions
diff --git a/src/templates.c b/src/templates.c index e4355295..943ddbd0 100644 --- a/src/templates.c +++ b/src/templates.c @@ -20,6 +20,9 @@ #include "utils.h" #include "geometry.h" #include "hdf5-file.h" +#include "peaks.h" + +#include <assert.h> /* Private data for template indexing */ @@ -35,7 +38,7 @@ struct template { double omega; double phi; int n; - struct reflhit spots; /* Made into an array by Magic */ + struct reflhit *spots; }; @@ -98,6 +101,7 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename, double omega, phi; struct image image; struct hdfile *hdfile; + int i; hdfile = hdfile_open(filename); if ( hdfile == NULL ) { @@ -135,13 +139,18 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename, n_templates = (omega_max * phi_max)/(omega_step * phi_step); STATUS("%i templates to be calculated.\n", n_templates); - for ( omega = 0.0; omega < omega_max; omega += omega_step ) { - for ( phi = 0.0; phi < phi_max; phi += phi_step ) { + priv->templates = malloc(n_templates * sizeof(struct template)); + + i = 0; + for ( omega = 0.0; omega < omega_max-omega_step; omega += omega_step ) { + for ( phi = 0.0; phi < phi_max-phi_step; phi += phi_step ) { int n; struct reflhit *hits; UnitCell *cell_rot; + assert(i < n_templates); + cell_rot = rotate_cell(cell, omega, phi); hits = find_intersections(&image, cell_rot, 5.0e-3, @@ -151,6 +160,12 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename, return NULL; } + priv->templates[i].omega = omega; + priv->templates[i].phi = phi; + priv->templates[i].n = n; + priv->templates[i].spots = hits; + i++; + free(cell_rot); } @@ -164,6 +179,66 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename, } +static double integrate_all_rot(struct image *image, struct reflhit *hits, + int n, double rot) +{ + double itot = 0.0; + int i; + + for ( i=0; i<n; i++ ) { + + float x, y, intensity; + float xp, yp; + + xp = cos(rot)*hits[i].x + sin(rot)*hits[i].y; + yp = -sin(rot)*hits[i].x + cos(rot)*hits[i].y; + + if ( integrate_peak(image, xp, yp, &x, &y, + &intensity, 0, 0) ) continue; + + itot += intensity; + } + + return itot; +} + + void match_templates(struct image *image, IndexingPrivate *ipriv) { + struct _indexingprivate_template *priv + = (struct _indexingprivate_template *)ipriv; + int i, max_i; + double max, tot; + double rot, rot_max, rot_step; + + max = 0.0; + max_i = 0; + rot_max = 2.0*M_PI; + rot_step = 2.0*M_PI / 360.0; /* 1 deg steps */ + + for ( i=0; i<priv->n_templates; i++ ) { + for ( rot=0.0; rot<rot_max; rot+=rot_step ) { + + double val; + + val = integrate_all_rot(image, priv->templates[i].spots, + priv->templates[i].n, rot); + + if ( val > max ) { + max = val; + max_i = i; + } + + } + progress_bar(i, priv->n_templates, "Indexing"); + } + + tot = 0.0; + for ( i=0; i<(image->width*image->height); i++ ) { + tot += image->data[i]; + } + + STATUS("%i (%.2f, %.2f): %.2f%%\n", max_i, priv->templates[max_i].omega, + priv->templates[max_i].phi, + 100.0 * max / tot); } |