From ac20d7cd290e7bac157b9b90a73a04a956ab6ee3 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Mon, 16 Aug 2010 17:10:11 +0200
Subject: Determine orientation range for template matching

---
 src/symmetry.c  | 87 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 src/symmetry.h  |  5 ++++
 src/templates.c | 25 +++++++++++++++++
 3 files changed, 117 insertions(+)

diff --git a/src/symmetry.c b/src/symmetry.c
index bb9e966e..6611e617 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -431,3 +431,90 @@ int find_unique_equiv(ReflItemList *items, signed int h, signed int k,
 
 	return found;
 }
+
+
+/* Returns true if the point group is 23, m-3, 432, -43m or m-3m
+ * (i.e. T, Th, O, Td or Oh) */
+int is_polyhedral(const char *sym)
+{
+	/* Triclinic */
+	if ( strcmp(sym, "1") == 0 ) return 0;
+	if ( strcmp(sym, "-1") == 0 ) return 0;
+
+	/* Tetragonal */
+	if ( strcmp(sym, "422") == 0 ) return 0;
+
+	/* Hexagonal */
+	if ( strcmp(sym, "6") == 0 ) return 0;
+	if ( strcmp(sym, "6/m") == 0 ) return 0;
+	if ( strcmp(sym, "6/mmm") == 0 ) return 0;
+
+	/* TODO: Add more groups here */
+
+	ERROR("Don't know if '%s' is polyhedral or not.\n", sym);
+	abort();
+}
+
+
+/* Returns the order of the highest axis of proper or improper rotation */
+int rotational_order(const char *sym)
+{
+	/* Triclinic */
+	if ( strcmp(sym, "1") == 0 ) return 1;
+	if ( strcmp(sym, "-1") == 0 ) return 2 ;
+
+	/* Tetragonal */
+	if ( strcmp(sym, "422") == 0 ) return 4;
+
+	/* Hexagonal */
+	if ( strcmp(sym, "6") == 0 ) return 6;
+	if ( strcmp(sym, "6/m") == 0 ) return 6;
+	if ( strcmp(sym, "6/mmm") == 0 ) return 6;
+
+	/* TODO: Add more groups here */
+
+	ERROR("Couldn't find rotational order for '%s'.\n", sym);
+	abort();
+}
+
+
+int has_perpendicular_mirror(const char *sym)
+{
+	/* Triclinic */
+	if ( strcmp(sym, "1") == 0 ) return 0;
+	if ( strcmp(sym, "-1") == 0 ) return 0;
+
+	/* Tetragonal */
+	if ( strcmp(sym, "422") == 0 ) return 0;
+
+	/* Hexagonal */
+	if ( strcmp(sym, "6") == 0 ) return 0;
+	if ( strcmp(sym, "6/m") == 0 ) return 1;
+	if ( strcmp(sym, "6/mmm") == 0 ) return 1;
+
+	/* TODO: Add more groups here */
+
+	ERROR("Couldn't find mirror definition for '%s'.\n", sym);
+	abort();
+}
+
+
+int has_bisecting_mirror_or_diad(const char *sym)
+{
+	/* Triclinic */
+	if ( strcmp(sym, "1") == 0 ) return 0;
+	if ( strcmp(sym, "-1") == 0 ) return 0;
+
+	/* Tetragonal */
+	if ( strcmp(sym, "422") == 0 ) return 0;
+
+	/* Hexagonal */
+	if ( strcmp(sym, "6") == 0 ) return 0;
+	if ( strcmp(sym, "6/m") == 0 ) return 1;
+	if ( strcmp(sym, "6/mmm") == 0 ) return 1;
+
+	/* TODO: Add more groups here */
+
+	ERROR("Couldn't find mirror definition for '%s'.\n", sym);
+	abort();
+}
diff --git a/src/symmetry.h b/src/symmetry.h
index bb4ffa60..c338da49 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -44,5 +44,10 @@ extern int find_unique_equiv(ReflItemList *items, signed int h, signed int k,
                              signed int l, const char *mero, signed int *hu,
                              signed int *ku, signed int *lu);
 
+/* Properties of point groups */
+extern int is_polyhedral(const char *sym);
+extern int rotational_order(const char *sym);
+extern int has_perpendicular_mirror(const char *sym);
+extern int has_bisecting_mirror_or_diad(const char *sym);
 
 #endif	/* SYMMETRY_H */
diff --git a/src/templates.c b/src/templates.c
index db28b9b9..98182bbe 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -16,6 +16,8 @@
 
 #include "index.h"
 #include "index-priv.h"
+#include "symmetry.h"
+#include "utils.h"
 
 
 /* Private data for template indexing */
@@ -25,13 +27,36 @@ struct _indexingprivate_template
 };
 
 
+/* Generate templates for the given cell using a representative image */
 IndexingPrivate *generate_templates(UnitCell *cell, const char *filename)
 {
 	struct _indexingprivate_template *priv;
+	const char *holo;
+	double omega_max, phi_max;
 
 	priv = calloc(1, sizeof(struct _indexingprivate_template));
 	priv->base.indm = INDEXING_TEMPLATE;
 
+	/* We can only distinguish orientations within the holohedral cell */
+	holo = get_holohedral(cell_get_pointgroup(cell));
+	STATUS("%s\n", holo);
+
+	/* These define the orientation in space */
+	if ( is_polyhedral(holo) ) {
+		ERROR("WARNING: Holohedral point group is polyhedral.\n");
+		ERROR("This means I can't properly determine the orientation");
+		ERROR(" ranges for template matching.  Expect trouble.\n");
+	}
+	omega_max = 2.0*M_PI / rotational_order(holo);
+	if ( has_bisecting_mirror_or_diad(holo) ) omega_max /= 2.0;
+	phi_max = M_PI;
+	if ( has_perpendicular_mirror(holo) ) phi_max /= 2.0;
+
+	/* One more axis would define the rotation in the plane of the image */
+
+	STATUS("Orientation ranges: %5.0f -> %5.0f, %5.0f -> %5.0f deg.\n",
+	       0.0, rad2deg(omega_max), 0.0, rad2deg(phi_max));
+
 	return (struct _indexingprivate *)priv;
 }
 
-- 
cgit v1.2.3