aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-08-17 18:14:44 +0200
committerThomas White <taw@physics.org>2012-02-22 15:26:55 +0100
commit0dcc9e2e1ea4fcbd37db5ab7ac74146098c2f4d7 (patch)
tree22fce4b990b75619c19613824289c3f71275635c
parentac20d7cd290e7bac157b9b90a73a04a956ab6ee3 (diff)
Generation of templates (needs debugging)
-rw-r--r--Makefile.am2
-rw-r--r--Makefile.in2
-rw-r--r--src/Makefile.am5
-rw-r--r--src/Makefile.in12
-rw-r--r--src/cell.c3
-rw-r--r--src/facetron.c135
-rw-r--r--src/symmetry.c3
-rw-r--r--src/templates.c109
8 files changed, 124 insertions, 147 deletions
diff --git a/Makefile.am b/Makefile.am
index e8cc6ffd..35c962e3 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -6,5 +6,5 @@ EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h \
src/filters.h src/diffraction-gpu.h src/cl-utils.h \
data/defs.h src/parameters-lcls.tmp \
data/diffraction.cl data/sfac src/likelihood.h src/symmetry.h \
- src/povray.h src/index-priv.h
+ src/povray.h src/index-priv.h src/geometry.h
SUBDIRS = src data
diff --git a/Makefile.in b/Makefile.in
index 55a85578..84406a2f 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -203,7 +203,7 @@ EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h \
src/filters.h src/diffraction-gpu.h src/cl-utils.h \
data/defs.h src/parameters-lcls.tmp \
data/diffraction.cl data/sfac src/likelihood.h src/symmetry.h \
- src/povray.h src/index-priv.h
+ src/povray.h src/index-priv.h src/geometry.h
SUBDIRS = src data
all: config.h
diff --git a/src/Makefile.am b/src/Makefile.am
index b7108918..7fdf6dc9 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -22,7 +22,8 @@ process_hkl_LDADD = @LIBS@
indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \
peaks.c index.c filters.c diffraction.c detector.c \
- sfac.c dirax.c reflections.c templates.c symmetry.c
+ sfac.c dirax.c reflections.c templates.c symmetry.c \
+ geometry.c
indexamajig_LDADD = @LIBS@
if HAVE_OPENCL
indexamajig_SOURCES += diffraction-gpu.c cl-utils.c
@@ -56,7 +57,7 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
calibrate_detector_LDADD = @LIBS@
facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
- image.c diffraction.c sfac.c
+ image.c diffraction.c sfac.c geometry.c
facetron_LDADD = @LIBS@
cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c
diff --git a/src/Makefile.in b/src/Makefile.in
index 9355ef39..4b4da1f7 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -71,7 +71,7 @@ cubeit_DEPENDENCIES =
am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \
hdf5-file.$(OBJEXT) utils.$(OBJEXT) detector.$(OBJEXT) \
peaks.$(OBJEXT) image.$(OBJEXT) diffraction.$(OBJEXT) \
- sfac.$(OBJEXT)
+ sfac.$(OBJEXT) geometry.$(OBJEXT)
facetron_OBJECTS = $(am_facetron_OBJECTS)
facetron_DEPENDENCIES =
am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \
@@ -89,7 +89,7 @@ hdfsee_DEPENDENCIES =
am__indexamajig_SOURCES_DIST = indexamajig.c hdf5-file.c utils.c \
cell.c image.c peaks.c index.c filters.c diffraction.c \
detector.c sfac.c dirax.c reflections.c templates.c symmetry.c \
- diffraction-gpu.c cl-utils.c
+ geometry.c diffraction-gpu.c cl-utils.c
@HAVE_OPENCL_TRUE@am__objects_1 = diffraction-gpu.$(OBJEXT) \
@HAVE_OPENCL_TRUE@ cl-utils.$(OBJEXT)
am_indexamajig_OBJECTS = indexamajig.$(OBJEXT) hdf5-file.$(OBJEXT) \
@@ -97,7 +97,7 @@ am_indexamajig_OBJECTS = indexamajig.$(OBJEXT) hdf5-file.$(OBJEXT) \
index.$(OBJEXT) filters.$(OBJEXT) diffraction.$(OBJEXT) \
detector.$(OBJEXT) sfac.$(OBJEXT) dirax.$(OBJEXT) \
reflections.$(OBJEXT) templates.$(OBJEXT) symmetry.$(OBJEXT) \
- $(am__objects_1)
+ geometry.$(OBJEXT) $(am__objects_1)
indexamajig_OBJECTS = $(am_indexamajig_OBJECTS)
indexamajig_DEPENDENCIES =
am__pattern_sim_SOURCES_DIST = pattern_sim.c diffraction.c utils.c \
@@ -260,7 +260,8 @@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \
process_hkl_LDADD = @LIBS@
indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \
peaks.c index.c filters.c diffraction.c detector.c sfac.c \
- dirax.c reflections.c templates.c symmetry.c $(am__append_3)
+ dirax.c reflections.c templates.c symmetry.c geometry.c \
+ $(am__append_3)
indexamajig_LDADD = @LIBS@
@HAVE_GTK_TRUE@hdfsee_SOURCES = hdfsee.c displaywindow.c render.c hdf5-file.c utils.c image.c \
@HAVE_GTK_TRUE@ filters.c
@@ -286,7 +287,7 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
calibrate_detector_LDADD = @LIBS@
facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
- image.c diffraction.c sfac.c
+ image.c diffraction.c sfac.c geometry.c
facetron_LDADD = @LIBS@
cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c
@@ -415,6 +416,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/displaywindow.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/facetron.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/filters.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geometry.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/get_hkl.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/hdf5-file.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/hdfsee.Po@am__quote@
diff --git a/src/cell.c b/src/cell.c
index d6464a71..d6d101cc 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -207,6 +207,9 @@ static int cell_crystallographic_to_cartesian(UnitCell *cell,
{
double tmp, V, cosalphastar, cstar;
+ /* Note: Please consider and possibly change the ranges for template
+ * matching (in templates.c) if the calculations below are altered. */
+
/* Firstly: Get a in terms of x, y and z
* +a (cryst) is defined to lie along +x (cart) */
*ax = cell->a;
diff --git a/src/facetron.c b/src/facetron.c
index 0fa6b643..0c2b4d56 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -26,10 +26,7 @@
#include "cell.h"
#include "hdf5-file.h"
#include "detector.h"
-#include "peaks.h"
-
-
-#define MAX_HITS (1024)
+#include "geometry.h"
static void show_help(const char *s)
@@ -116,136 +113,6 @@ static int find_chunk(FILE *fh, UnitCell **cell, char **filename)
}
-static struct reflhit *find_intersections(struct image *image, UnitCell *cell,
- double divergence, double bandwidth,
- int *n, int output)
-{
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- struct reflhit *hits;
- int np = 0;
- int hmax, kmax, lmax;
- double mres;
- signed int h, k, l;
-
- hits = malloc(sizeof(struct reflhit)*MAX_HITS);
- if ( hits == NULL ) {
- *n = 0;
- return NULL;
- }
-
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
-
- mres = 1.0 / 8.0e-10; /* 8 Angstroms */
- hmax = mres / modulus(asx, asy, asz);
- kmax = mres / modulus(bsx, bsy, bsz);
- lmax = mres / modulus(csx, csy, csz);
-
- for ( h=-hmax; h<hmax; h++ ) {
- for ( k=-kmax; k<kmax; k++ ) {
- for ( l=-lmax; l<lmax; l++ ) {
-
- double xl, yl, zl;
- double ds_sq, dps_sq;
- double delta, divfact;
- double llow, lhigh;
- double xd, yd, cl;
- double xda, yda;
- int p;
- int found = 0;
-
- if ( (h==0) && (k==0) && (l==0) ) continue;
-
- llow = image->lambda - image->lambda*bandwidth/2.0;
- lhigh = image->lambda + image->lambda*bandwidth/2.0;
-
- /* Get the coordinates of the reciprocal lattice point */
- zl = h*asz + k*bsz + l*csz;
- if ( zl < 0.0 ) continue; /* Do this check very early */
- xl = h*asx + k*bsx + l*csx;
- yl = h*asy + k*bsy + l*csy;
-
- ds_sq = modulus_squared(xl, yl, zl); /* d*^2 */
- delta = divergence/image->lambda;
- dps_sq = ds_sq + pow(delta, 2.0); /* d'*^2 */
-
- /* In range? */
- divfact = 2.0 * delta * sqrt(xl*xl + yl*yl);
- if ( ds_sq - 2.0*zl/llow > 0.0 ) continue;
- if ( ds_sq - 2.0*zl/lhigh < 0.0 ) continue;
-
- /* Work out which panel this peak would fall on */
- for ( p=0; p<image->det->n_panels; p++ ) {
-
- /* Camera length for this panel */
- cl = image->det->panels[p].clen;
-
- /* Coordinates of peak relative to central beam, in m */
- xd = cl*xl / (ds_sq/(2.0*zl) - zl);
- yd = cl*yl / (ds_sq/(2.0*zl) - zl);
-
- /* Convert to pixels */
- xd *= image->det->panels[p].res;
- yd *= image->det->panels[p].res;
-
- /* Add the coordinates of the central beam */
- xda = xd + image->det->panels[p].cx;
- yda = yd + image->det->panels[p].cy;
-
- /* Now, is this on this panel? */
- if ( xda < image->det->panels[p].min_x ) continue;
- if ( xda > image->det->panels[p].max_x ) continue;
- if ( yda < image->det->panels[p].min_y ) continue;
- if ( yda > image->det->panels[p].max_y ) continue;
-
- /* Woohoo! */
- found = 1;
- break;
-
- }
-
- if ( !found ) continue;
-
- hits[np].h = h;
- hits[np].k = k;
- hits[np].l = l;
- np++;
-
- if ( output ) {
- printf("%i %i %i 0.0 (at %f,%f)\n", h, k, l, xda, yda);
- }
-
- }
- }
- }
-
- *n = np;
- return hits;
-}
-
-
-static double integrate_all(struct image *image, struct reflhit *hits, int n)
-{
- double itot = 0.0;
- int i;
-
- for ( i=0; i<n; i++ ) {
-
- float x, y, intensity;
-
- if ( integrate_peak(image, hits[i].x, hits[i].y, &x, &y,
- &intensity, 0, 0) ) continue;
-
- itot += intensity;
- }
-
- return itot;
-}
-
-
static void get_trial_cell(UnitCell *out, UnitCell *in, int i, double step)
{
double asx, asy, asz;
diff --git a/src/symmetry.c b/src/symmetry.c
index 6611e617..6511d173 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -456,7 +456,8 @@ int is_polyhedral(const char *sym)
}
-/* Returns the order of the highest axis of proper or improper rotation */
+/* Returns the order of the characteristic axis of proper or improper rotation
+ * for the point group. This should be the rotation about the c axis. */
int rotational_order(const char *sym)
{
/* Triclinic */
diff --git a/src/templates.c b/src/templates.c
index 98182bbe..fd4616b2 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -18,28 +18,101 @@
#include "index-priv.h"
#include "symmetry.h"
#include "utils.h"
+#include "geometry.h"
+#include "hdf5-file.h"
/* Private data for template indexing */
struct _indexingprivate_template
{
struct _indexingprivate base;
+ int n_templates;
+ struct template *templates;
};
+struct template {
+ double omega;
+ double phi;
+ int n;
+ struct reflhit spots; /* Made into an array by Magic */
+};
+
+
+UnitCell *rotate_cell(UnitCell *in, double omega, double phi)
+{
+ UnitCell *out;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double xnew, ynew, znew;
+
+ cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy,
+ &bsz, &csx, &csy, &csz);
+
+ /* Rotate by "omega" about +z (parallel to c* and c unless triclinic) */
+ xnew = asx*cos(omega) + asy*sin(omega);
+ ynew = -asx*sin(omega) + asy*cos(omega);
+ znew = asz;
+ asx = xnew; asy = ynew; asz = znew;
+ xnew = bsx*cos(omega) + bsy*sin(omega);
+ ynew = -bsx*sin(omega) + bsy*cos(omega);
+ znew = bsz;
+ bsx = xnew; bsy = ynew; bsz = znew;
+ xnew = csx*cos(omega) + csy*sin(omega);
+ ynew = -csx*sin(omega) + csy*cos(omega);
+ znew = csz;
+ csx = xnew; csy = ynew; csz = znew;
+
+ /* Rotate by "phi" about +x (not parallel to anything specific) */
+ xnew = asx;
+ ynew = asy*cos(phi) + asz*sin(phi);
+ znew = -asy*sin(phi) + asz*cos(phi);
+ asx = xnew; asy = ynew; asz = znew;
+ xnew = bsx;
+ ynew = bsy*cos(phi) + bsz*sin(phi);
+ znew = -bsy*sin(phi) + bsz*cos(phi);
+ bsx = xnew; bsy = ynew; bsz = znew;
+ xnew = csx;
+ ynew = csy*cos(phi) + csz*sin(phi);
+ znew = -csy*sin(phi) + csz*cos(phi);
+ csx = xnew; csy = ynew; csz = znew;
+
+ out = cell_new_from_cell(in);
+ cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+
+ return out;
+}
+
+
/* 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;
+ int n_templates;
+ const double omega_step = deg2rad(0.5);
+ const double phi_step = deg2rad(0.5);
+ double omega, phi;
+ struct image image;
+ struct hdfile *hdfile;
+
+ hdfile = hdfile_open(filename);
+ if ( hdfile == NULL ) {
+ return NULL;
+ } else if ( hdfile_set_image(hdfile, "/data/data0") ) {
+ ERROR("Couldn't select path\n");
+ return NULL;
+ }
+ hdf5_read(hdfile, &image, 0);
+ hdfile_close(hdfile);
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) ) {
@@ -54,8 +127,38 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename)
/* 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));
+ STATUS("Orientation ranges in %s: %.0f-%.0f, %.0f-%.0f deg.\n",
+ holo, 0.0, rad2deg(omega_max), 0.0, rad2deg(phi_max));
+
+ 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 ) {
+
+ int n;
+ struct reflhit *hits;
+ UnitCell *cell_rot;
+
+ cell_rot = rotate_cell(cell, omega, phi);
+
+ hits = find_intersections(&image, cell_rot, 5.0e-3,
+ 3.0/100.0, &n, 0);
+ if ( hits == NULL ) {
+ ERROR("Template calculation failed.\n");
+ return NULL;
+ }
+
+ STATUS("%.2f, %.2f : %i features\n",
+ rad2deg(omega), rad2deg(phi), n);
+
+ free(cell_rot);
+
+ }
+ STATUS("Finished omega=%.2f\n", rad2deg(omega));
+ }
+
+ priv->n_templates = n_templates;
return (struct _indexingprivate *)priv;
}