aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-08-20 15:45:18 +0200
committerThomas White <taw@physics.org>2012-02-22 15:26:56 +0100
commit2e563d59a7f3f93b4fe1b1170dc6b47043c7170c (patch)
tree7b4c5f4f38b1bd63009cf8241c7e6e59376cc030 /src
parent13310cf31ef3aa731ab76fc3f7880a6da085e50e (diff)
Implement template indexing
Slllloooowwww......!
Diffstat (limited to 'src')
-rw-r--r--src/templates.c81
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);
}