Image, feature and unit cell infrastructure

Brought across from DTR and Synth3D

12 files changed, 959 insertions, 2 deletions

diff --git a/.gitignore b/.gitignore
index 3c13d553..692c0da4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,3 +10,4 @@ stamp-h1
 src/*.o
 src/.deps
 src/template_index
+src/simulate_patterns
diff --git a/src/Makefile.am b/src/Makefile.am
index ac173bb1..b480cedd 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,8 @@
-bin_PROGRAMS = template_index
-template_index_SOURCES = main.c
+bin_PROGRAMS = template_index simulate_patterns
+
+template_index_SOURCES = main.c relrod.c utils.c image.c cell.c
 template_index_LDADD = @LIBS@ -lm
 AM_CFLAGS = -Wall -g
+
+simulate_patterns_SOURCES = sim-main.c relrod.c utils.c image.c cell.c
+simulate_patterns_LDADD = @LIBS@ -lm
diff --git a/src/cell.c b/src/cell.c
new file mode 100644
index 00000000..d2495b0f
--- /dev/null
+++ b/src/cell.c
@@ -0,0 +1,148 @@
+/*
+ * cell.c
+ *
+ * Unit Cell Calculations
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "cell.h"
+#include "utils.h"
+
+
+/* Update the cartesian representation from the crystallographic one */
+static void cell_update_cartesian(UnitCell *cell)
+{
+	double tmp, V, cosalphastar, cstar;
+
+	if ( !cell ) return;
+
+	/* a in terms of x, y and z
+	 * +a (cryst) is defined to lie along +x (cart) */
+	cell->ax = cell->a;
+	cell->ay = 0.0;
+	cell->az = 0.0;
+
+	/* b in terms of x, y and z
+	 * b (cryst) is defined to lie in the xy (cart) plane */
+	cell->bx = cell->b*cos(cell->gamma);
+	cell->by = cell->b*sin(cell->gamma);
+	cell->bz = 0.0;
+
+	tmp = cos(cell->alpha)*cos(cell->alpha)
+		- cos(cell->beta)*cos(cell->beta)
+		- cos(cell->gamma)*cos(cell->gamma)
+		+ 2.0*cos(cell->alpha)*cos(cell->beta)*cos(cell->gamma);
+	V = cell->a * cell->b * cell->c * sqrt(1.0 - tmp);
+
+	cosalphastar = cos(cell->beta)*cos(cell->gamma) - cos(cell->alpha);
+	cosalphastar /= sin(cell->beta)*sin(cell->gamma);
+
+	cstar = (cell->a * cell->b * sin(cell->gamma))/V;
+
+	/* c in terms of x, y and z */
+	cell->cx = cell->c*cos(cell->beta);
+	cell->cy = -cell->c*sin(cell->beta)*cosalphastar;
+	cell->cz = 1.0/cstar;
+}
+
+
+/* Update the crystallographic representation from the cartesian one */
+static void cell_update_crystallographic(UnitCell *cell)
+{
+	if ( !cell ) return;
+
+	cell->a = modulus(cell->ax, cell->ay, cell->az);
+	cell->b = modulus(cell->bx, cell->by, cell->bz);
+	cell->c = modulus(cell->cx, cell->cy, cell->cz);
+
+	cell->alpha = angle_between(cell->bx, cell->by, cell->bz,
+					cell->cx, cell->cy, cell->cz);
+	cell->beta = angle_between(cell->ax, cell->ay, cell->az,
+					cell->cx, cell->cy, cell->cz);
+	cell->gamma = angle_between(cell->ax, cell->ay, cell->az,
+					cell->bx, cell->by, cell->bz);
+
+	printf("a=%f nm\n", cell->a);
+	printf("b=%f nm\n", cell->b);
+	printf("c=%f nm\n", cell->c);
+	printf("alpha = %f deg\n", rad2deg(cell->alpha));
+	printf(" beta = %f deg\n", rad2deg(cell->beta));
+	printf("gamma = %f deg\n", rad2deg(cell->gamma));
+}
+
+
+UnitCell *cell_new()
+{
+	UnitCell *cell;
+
+	cell = malloc(sizeof(UnitCell));
+	if ( !cell ) return NULL;
+
+	cell->a = 1.0;
+	cell->b = 1.0;
+	cell->c = 1.0;
+	cell->alpha = M_PI_2;
+	cell->beta = M_PI_2;
+	cell->gamma = M_PI_2;
+
+	cell_update_cartesian(cell);
+
+	return cell;
+}
+
+
+void cell_set_parameters(UnitCell *cell, double a, double b, double c,
+				double alpha, double beta, double gamma)
+{
+	if ( !cell ) return;
+
+	cell->a = a;
+	cell->b = b;
+	cell->c = c;
+	cell->alpha = alpha;
+	cell->beta = beta;
+	cell->gamma = gamma;
+
+	cell_update_cartesian(cell);
+}
+
+
+void cell_set_cartesian(UnitCell *cell,
+			double ax, double ay, double az,
+			double bx, double by, double bz,
+			double cx, double cy, double cz)
+{
+	if ( !cell ) return;
+
+	cell->ax = ax;  cell->ay = ay;  cell->az = az;
+	cell->bx = bx;  cell->by = by;  cell->bz = bz;
+	cell->cx = cx;  cell->cy = cy;  cell->cz = cz;
+
+	cell_update_crystallographic(cell);
+}
+
+
+UnitCell *cell_new_from_parameters(double a, double b, double c,
+				double alpha, double beta, double gamma)
+{
+	UnitCell *cell;
+
+	cell = cell_new();
+	if ( !cell ) return NULL;
+
+	cell_set_parameters(cell, a, b, c, alpha, beta, gamma);
+
+	return cell;
+}
diff --git a/src/cell.h b/src/cell.h
new file mode 100644
index 00000000..e99c88df
--- /dev/null
+++ b/src/cell.h
@@ -0,0 +1,55 @@
+/*
+ * cell.h
+ *
+ * Unit Cell Calculations
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+#ifndef CELL_H
+#define CELL_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+typedef struct {
+
+	/* Crystallographic representation */
+	double a;	/* nm */
+	double b;	/* nm */
+	double c;	/* nm */
+	double alpha;	/* Radians */
+	double beta;	/* Radians */
+	double gamma;	/* Radians */
+
+	/* Cartesian representation */
+	double ax;	double bx;	double cx;
+	double ay;	double by;	double cy;
+	double az;	double bz;	double cz;
+
+	/* Cartesian representation of reciprocal axes */
+	double axs;	double bxs;	double cxs;
+	double ays;	double bys;	double cys;
+	double azs;	double bzs;	double czs;
+
+} UnitCell;
+
+extern UnitCell *cell_new(void);
+
+/* Lengths in nm, angles in radians */
+extern UnitCell *cell_new_from_parameters(double a, double b, double c,
+				double alpha, double beta, double gamma);
+
+extern void cell_set_cartesian(UnitCell *cell,
+			double ax, double ay, double az,
+			double bx, double by, double bz,
+			double cx, double cy, double cz);
+
+extern void cell_set_parameters(UnitCell *cell, double a, double b, double c,
+				double alpha, double beta, double gamma);
+
+#endif	/* CELL_H */
diff --git a/src/image.c b/src/image.c
new file mode 100644
index 00000000..b3825b38
--- /dev/null
+++ b/src/image.c
@@ -0,0 +1,174 @@
+/*
+ * image.c
+ *
+ * Handle images and image features
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+
+#define _GNU_SOURCE 1
+#include <stdlib.h>
+#include <assert.h>
+#include <math.h>
+
+#include "image.h"
+#include "utils.h"
+
+
+int image_add(ImageList *list, struct image *image)
+{
+	if ( list->images ) {
+		list->images = realloc(list->images,
+		                       (list->n_images+1)*sizeof(struct image));
+	} else {
+		assert(list->n_images == 0);
+		list->images = malloc(sizeof(struct image));
+	}
+
+	/* Copy the metadata */
+	list->images[list->n_images] = *image;
+
+	/* Fill out some extra fields */
+	list->images[list->n_images].features = NULL;
+	list->images[list->n_images].rflist = NULL;
+
+	list->n_images++;
+
+	return list->n_images - 1;
+}
+
+
+ImageList *image_list_new()
+{
+	ImageList *list;
+
+	list = malloc(sizeof(ImageList));
+
+	list->n_images = 0;
+	list->images = NULL;
+
+	return list;
+}
+
+
+void image_add_feature_reflection(ImageFeatureList *flist, double x, double y,
+                                  struct image *parent, double intensity)
+{
+	if ( flist->features ) {
+		flist->features = realloc(flist->features,
+		                    (flist->n_features+1)*sizeof(ImageFeature));
+	} else {
+		assert(flist->n_features == 0);
+		flist->features = malloc(sizeof(ImageFeature));
+	}
+
+	flist->features[flist->n_features].x = x;
+	flist->features[flist->n_features].y = y;
+	flist->features[flist->n_features].intensity = intensity;
+	flist->features[flist->n_features].parent = parent;
+	flist->features[flist->n_features].partner = NULL;
+	flist->features[flist->n_features].partner_d = 0.0;
+
+	flist->n_features++;
+
+}
+
+
+void image_add_feature(ImageFeatureList *flist, double x, double y,
+                       struct image *parent, double intensity)
+{
+	image_add_feature_reflection(flist, x, y, parent, intensity);
+}
+
+
+ImageFeatureList *image_feature_list_new()
+{
+	ImageFeatureList *flist;
+
+	flist = malloc(sizeof(ImageFeatureList));
+
+	flist->n_features = 0;
+	flist->features = NULL;
+
+	return flist;
+}
+
+
+void image_feature_list_free(ImageFeatureList *flist)
+{
+	if ( !flist ) return;
+
+	if ( flist->features ) free(flist->features);
+	free(flist);
+}
+
+
+ImageFeature *image_feature_closest(ImageFeatureList *flist, double x, double y,
+                                    double *d, int *idx)
+{
+	int i;
+	double dmin = +HUGE_VAL;
+	int closest = 0;
+
+	for ( i=0; i<flist->n_features; i++ ) {
+
+		double d;
+
+		d = distance(flist->features[i].x, flist->features[i].y, x, y);
+
+		if ( d < dmin ) {
+			dmin = d;
+			closest = i;
+		}
+
+	}
+
+	if ( dmin < +HUGE_VAL ) {
+		*d = dmin;
+		*idx = closest;
+		return &flist->features[closest];
+	}
+
+	*d = +INFINITY;
+	return NULL;
+}
+
+
+ImageFeature *image_feature_second_closest(ImageFeatureList *flist,
+                                           double x, double y, double *d,
+                                           int *idx)
+{
+	int i;
+	double dmin = +HUGE_VAL;
+	int closest = 0;
+	double dfirst;
+	int idxfirst;
+
+	image_feature_closest(flist, x, y, &dfirst, &idxfirst);
+
+	for ( i=0; i<flist->n_features; i++ ) {
+
+		double d;
+
+		d = distance(flist->features[i].x, flist->features[i].y, x, y);
+
+		if ( (d < dmin) && (i != idxfirst) ) {
+			dmin = d;
+			closest = i;
+		}
+
+	}
+
+	if ( dmin < +HUGE_VAL ) {
+		*d = dmin;
+		*idx = closest;
+		return &flist->features[closest];
+	}
+
+	return NULL;
+
+}
diff --git a/src/image.h b/src/image.h
new file mode 100644
index 00000000..f53f7d92
--- /dev/null
+++ b/src/image.h
@@ -0,0 +1,121 @@
+/*
+ * image.h
+ *
+ * Handle images and image features
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef IMAGE_H
+#define IMAGE_H
+
+#include <stdint.h>
+
+
+typedef enum {
+	FORMULATION_CLEN,
+	FORMULATION_PIXELSIZE
+} FormulationMode;
+
+
+typedef struct imagefeature_struct {
+
+	struct image			*parent;
+	double				x;
+	double				y;
+	double				intensity;
+
+	/* Partner for this feature (in another feature list) or NULL */
+	struct imagefeature_struct	*partner;
+
+	/* Distance between this feature and its partner, if any. */
+	double				partner_d;
+
+	/* The reflection this was projected from, if any */
+	struct reflection_struct	*reflection;
+
+} ImageFeature;
+
+
+typedef struct {
+
+	ImageFeature		*features;
+	int			n_features;
+
+} ImageFeatureList;
+
+
+struct image {
+
+	uint16_t		*data;
+
+	/* Radians.  Defines where the pattern lies in reciprocal space */
+	double			tilt;
+
+	/* Radians.  Defines where the pattern lies in reciprocal space */
+	double			omega;
+
+	/* Image parameters can be given as camera length or pixel size.
+	 * If FORMULATION_CLEN, then camera_len and resolution must be given.
+	 * If FORMULATION_PIXELSIZE, then pixel_size only is needed.*/
+	FormulationMode		fmode;
+	double			pixel_size;
+	double			camera_len;
+	double			resolution;	/* pixels per metre */
+
+	/* Wavelength must always be given */
+	double			lambda;
+
+	int			width;
+	int			height;
+	double			x_centre;
+	double			y_centre;
+
+	ImageFeatureList	*features;	/* "Experimental" features */
+	ImageFeatureList	*rflist;	/* "Predicted" features */
+
+};
+
+
+typedef struct imagelist_struct {
+
+	int		n_images;
+	struct image	*images;
+
+} ImageList;
+
+
+extern ImageList *image_list_new(void);
+
+extern int image_add(ImageList *list, struct image *image);
+
+extern ImageFeatureList *image_feature_list_new(void);
+
+extern void image_feature_list_free(ImageFeatureList *flist);
+
+extern void image_add_feature(ImageFeatureList *flist, double x, double y,
+                              struct image *parent, double intensity);
+
+extern void image_add_feature_reflection(ImageFeatureList *flist,
+                                         double x, double y,
+                                         struct image *parent,
+                                         double intensity);
+
+extern ImageFeature *image_feature_closest(ImageFeatureList *flist,
+                                           double x, double y, double *d,
+                                           int *idx);
+
+extern ImageFeature *image_feature_second_closest(ImageFeatureList *flist,
+                                                  double x, double y, double *d,
+                                                  int *idx);
+
+
+#endif	/* IMAGE_H */
diff --git a/src/main.c b/src/main.c
index ffb269e4..79dff45a 100644
--- a/src/main.c
+++ b/src/main.c
@@ -60,6 +60,8 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	printf("Generating templates...\n");
+
 	printf("Input files (%i):\n", nin);
 	for ( i=0; i<nin; i++ ) {
 		printf("%6i: %s\n", i+1, in_files[i]);
diff --git a/src/relrod.c b/src/relrod.c
new file mode 100644
index 00000000..674f7971
--- /dev/null
+++ b/src/relrod.c
@@ -0,0 +1,180 @@
+/*
+ * relrod.c
+ *
+ * Calculate reflection positions via line-sphere intersection test
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "image.h"
+#include "utils.h"
+#include "cell.h"
+
+
+static void mapping_rotate(double x, double y, double z,
+                           double *ddx, double *ddy, double *ddz,
+                           double omega, double tilt)
+{
+	double nx, ny, nz;
+	double x_temp, y_temp, z_temp;
+
+	/* First: rotate image clockwise until tilt axis is aligned
+	 * horizontally. */
+	nx = x*cos(omega) + y*sin(omega);
+	ny = -x*sin(omega) + y*cos(omega);
+	nz = z;
+
+	/* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the
+	 * "wrong" way. This is because the crystal is rotated in the
+	 * experiment, not the Ewald sphere. */
+	x_temp = nx; y_temp = ny; z_temp = nz;
+	nx = x_temp;
+	ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
+	nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
+
+	/* Finally, reverse the omega rotation to restore the location of the
+	 * image in 3D space */
+	x_temp = nx; y_temp = ny; z_temp = nz;
+	nx = x_temp*cos(-omega) + y_temp*sin(-omega);
+	ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
+	nz = z_temp;
+
+	*ddx = nx;
+	*ddy = ny;
+	*ddz = nz;
+}
+
+
+void get_reflections(struct image *image, UnitCell *cell)
+{
+	ImageFeatureList *flist;
+	double smax = 0.2e9;
+	double tilt, omega, k;
+	double nx, ny, nz; /* "normal" vector */
+	double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda */
+	double ux, uy, uz; /* "up" vector */
+	double rx, ry, rz; /* "right" vector */
+
+	flist = image_feature_list_new();
+
+	tilt = image->tilt;
+	omega = image->omega;
+	k = 1.0/image->lambda;
+
+	/* Calculate the (normalised) incident electron wavevector */
+	mapping_rotate(0.0, 0.0, 1.0, &nx, &ny, &nz, omega, tilt);
+	kx = nx / image->lambda;
+	ky = ny / image->lambda;
+	kz = nz / image->lambda;	/* This is the centre of the Ewald sphere */
+
+	/* Determine where "up" is */
+	mapping_rotate(0.0, 1.0, 0.0, &ux, &uy, &uz, omega, tilt);
+
+	/* Determine where "right" is */
+	mapping_rotate(1.0, 0.0, 0.0, &rx, &ry, &rz, omega, tilt);
+
+	do {
+
+		double xl, yl, zl;
+		double a, b, c;
+		double s1, s2, s, t;
+		double g_sq, gn;
+
+		/* Get the coordinates of the reciprocal lattice point */
+	//	xl = reflection->x;
+	//	yl = reflection->y;
+	//	zl = reflection->z;
+		g_sq = modulus_squared(xl, yl, zl);
+		gn = xl*nx + yl*ny + zl*nz;
+
+		/* Next, solve the relrod equation to calculate
+		 * the excitation error */
+		a = 1.0;
+		b = 2.0*(gn - k);
+		c = -2.0*gn*k + g_sq;
+		t = -0.5*(b + sign(b)*sqrt(b*b - 4.0*a*c));
+		s1 = t/a;
+		s2 = c/t;
+		if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2;
+
+		/* Skip this reflection if s is large */
+		if ( fabs(s) <= smax ) {
+
+			double xi, yi, zi;
+			double gx, gy, gz;
+			double theta;
+			double x, y;
+
+			/* Determine the intersection point */
+			xi = xl + s*nx;  yi = yl + s*ny;  zi = zl + s*nz;
+
+			/* Calculate Bragg angle */
+			gx = xi - kx;
+			gy = yi - ky;
+			gz = zi - kz;	/* This is the vector from the centre of
+					*  the sphere to the intersection */
+			theta = angle_between(-kx, -ky, -kz, gx, gy, gz);
+
+			if ( theta > 0.0 ) {
+
+				double dx, dy, psi;
+
+				/* Calculate azimuth of point in image
+				 * (anticlockwise from +x) */
+				dx = xi*rx + yi*ry + zi*rz;
+				dy = xi*ux + yi*uy + zi*uz;
+				psi = atan2(dy, dx);
+
+				/* Get image coordinates from polar
+				 * representation */
+				if ( image->fmode == FORMULATION_CLEN ) {
+					x = image->camera_len*tan(theta)*cos(psi);
+					y = image->camera_len*tan(theta)*sin(psi);
+					x *= image->resolution;
+					y *= image->resolution;
+				} else if ( image->fmode==FORMULATION_PIXELSIZE ) {
+					x = tan(theta)*cos(psi) / image->lambda;
+					y = tan(theta)*sin(psi) / image->lambda;
+					x /= image->pixel_size;
+					y /= image->pixel_size;
+				} else {
+					fprintf(stderr,
+						"Unrecognised formulation mode "
+						"in reproject_get_reflections\n");
+					return NULL;
+				}
+
+				x += image->x_centre;
+				y += image->y_centre;
+
+				/* Sanity check */
+				if ( (x>=0) && (x<image->width)
+				  && (y>=0) && (y<image->height) ) {
+
+					/* Record the reflection
+					 *  Intensity should be multiplied by
+					 *  relrod spike function except
+					 *  reprojected reflections aren't used
+					 *  quantitatively for anything
+					 */
+					image_add_feature_reflection(flist, x,
+						y, image, 1.0);
+
+				} /* else it's outside the picture somewhere */
+
+			} /* else this is the central beam */
+
+		}
+
+	} while ( 1 );
+
+	image->rflist = flist;
+}
diff --git a/src/relrod.h b/src/relrod.h
new file mode 100644
index 00000000..1dac6057
--- /dev/null
+++ b/src/relrod.h
@@ -0,0 +1,24 @@
+/*
+ * relrod.h
+ *
+ * Calculate reflection positions via line-sphere intersection test
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef RELROD_H
+#define RELROD_H
+
+#include "image.h"
+#include "cell.h"
+
+extern void get_reflections(struct image *image, UnitCell *cell);
+
+#endif	/* RELROD_H */
diff --git a/src/sim-main.c b/src/sim-main.c
new file mode 100644
index 00000000..38ad45bc
--- /dev/null
+++ b/src/sim-main.c
@@ -0,0 +1,82 @@
+/*
+ * sim-main.c
+ *
+ * Simulate test data
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+
+#include "image.h"
+#include "relrod.h"
+#include "cell.h"
+#include "utils.h"
+
+
+static void main_show_help(const char *s)
+{
+	printf("Syntax: %s <file1.h5> <file2.h5> [...]\n\n", s);
+	printf("Index diffraction patterns\n\n");
+	printf("  -h              Display this help message\n");
+}
+
+
+int main(int argc, char *argv[])
+{
+	int c;
+	ImageList *list;
+	UnitCell *cell;
+	struct image image;
+
+	while ((c = getopt(argc, argv, "h")) != -1) {
+
+		switch ( c ) {
+
+			case 'h' : {
+				main_show_help(argv[0]);
+				return 0;
+			}
+
+		}
+
+	}
+
+	printf("Generating test data...\n");
+	list = image_list_new();
+	image.width = 512;
+	image.height = 512;
+	image.tilt = 0.0;
+	image.omega = 0.0;
+	image.fmode = FORMULATION_CLEN;
+	image.x_centre = 128;
+	image.y_centre = 128;
+	image.camera_len = 1.0; /* one metre */
+	image.resolution = 5120;
+	image.lambda = 0.6e-9;
+	image.data = malloc(512*512*2);
+	image_add(list, &image);
+
+	cell = cell_new_from_parameters(1.0,
+	                                1.0,
+	                                1.0,
+	                                deg2rad(90.0),
+	                                deg2rad(90.0),
+	                                deg2rad(90.0));
+
+	get_reflections(&image, cell);
+
+	return 0;
+}
diff --git a/src/utils.c b/src/utils.c
new file mode 100644
index 00000000..0c7d7ad4
--- /dev/null
+++ b/src/utils.c
@@ -0,0 +1,124 @@
+/*
+ * utils.c
+ *
+ * Utility stuff
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+#include <math.h>
+#include <string.h>
+
+#include "utils.h"
+
+
+/* Return the MOST POSITIVE of two numbers */
+unsigned int biggest(signed int a, signed int b)
+{
+	if ( a>b ) {
+		return a;
+	}
+	return b;
+}
+
+
+/* Return the LEAST POSITIVE of two numbers */
+unsigned int smallest(signed int a, signed int b)
+{
+	if ( a<b ) {
+		return a;
+	}
+	return b;
+}
+
+
+double distance(double x1, double y1, double x2, double y2)
+{
+	return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
+}
+
+
+double modulus(double x, double y, double z)
+{
+	return sqrt(x*x + y*y + z*z);
+}
+
+
+double modulus_squared(double x, double y, double z) {
+	return x*x + y*y + z*z;
+}
+
+
+double distance3d(double x1, double y1, double z1,
+                  double x2, double y2, double z2)
+{
+	return modulus(x1-x2, y1-y2, z1-z2);
+}
+
+
+/* Angle between two vectors.  Answer in radians */
+double angle_between(double x1, double y1, double z1,
+                     double x2, double y2, double z2)
+{
+	double mod1 = modulus(x1, y1, z1);
+	double mod2 = modulus(x2, y2, z2);
+	return acos( (x1*x2 + y1*y2 + z1*z2) / (mod1*mod2) );
+}
+
+
+/* As above, answer in degrees */
+double angle_between_d(double x1, double y1, double z1,
+                       double x2, double y2, double z2)
+{
+	return rad2deg(angle_between(x1, y1, z1, x2, y2, z2));
+}
+
+
+/* Wavelength of an electron (in m) given accelerating potential (in V) */
+double lambda(double V)
+{
+	double m = 9.110E-31;
+	double h = 6.625E-34;
+	double e = 1.60E-19;
+	double c = 2.998E8;
+
+	return h / sqrt(2*m*e*V*(1+((e*V) / (2*m*c*c))));
+}
+
+
+size_t skipspace(const char *s)
+{
+	size_t i;
+
+	for ( i=0; i<strlen(s); i++ ) {
+		if ( (s[i] != ' ') && (s[i] != '\t') ) return i;
+	}
+
+	return strlen(s);
+}
+
+
+void chomp(char *s)
+{
+	size_t i;
+
+	if ( !s ) return;
+
+	for ( i=0; i<strlen(s); i++ ) {
+		if ( (s[i] == '\n') || (s[i] == '\r') ) {
+			s[i] = '\0';
+			return;
+		}
+	}
+}
+
+
+int sign(double a)
+{
+	if ( a < 0 ) return -1;
+	if ( a > 0 ) return +1;
+	return 0;
+}
diff --git a/src/utils.h b/src/utils.h
new file mode 100644
index 00000000..ef37bc73
--- /dev/null
+++ b/src/utils.h
@@ -0,0 +1,42 @@
+/*
+ * utils.h
+ *
+ * Utility stuff
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+#ifndef UTILS_H
+#define UTILS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+
+extern unsigned int biggest(signed int a, signed int b);
+extern unsigned int smallest(signed int a, signed int b);
+extern double distance(double x1, double y1, double x2, double y2);
+extern double modulus(double x, double y, double z);
+extern double modulus_squared(double x, double y, double z);
+extern double angle_between(double x1, double y1, double z1,
+                            double x2, double y2, double z2);
+extern double angle_between_d(double x1, double y1, double z1,
+                              double x2, double y2, double z2);
+extern double lambda(double voltage);
+extern double distance3d(double x1, double y1, double z1,
+                         double x2, double y2, double z2);
+extern size_t skipspace(const char *s);
+extern void chomp(char *s);
+extern int sign(double a);
+
+#define rad2deg(a) ((a)*180/M_PI)
+#define deg2rad(a) ((a)*M_PI/180)
+
+#define is_odd(a) ((a)%2==1)
+
+#endif	/* UTILS_H */