Read image and generate templates

Makefile.am

-EXTRA_DIST = configure src/reflist.h src/statistics.h
+EXTRA_DIST = configure src/cell.h src/hdf5-file.h src/image.h src/relrod.h \
+             src/templates.h src/utils.h
 SUBDIRS = src

src/Makefile.am

 bin_PROGRAMS = template_index simulate_patterns
 
-template_index_SOURCES = main.c relrod.c utils.c image.c cell.c hdf5-file.c
+template_index_SOURCES = main.c relrod.c utils.c image.c cell.c hdf5-file.c \
+                         templates.c
 template_index_LDADD = @LIBS@ -lm -lgsl -lgslcblas
 AM_CFLAGS = -Wall -g
 
 simulate_patterns_SOURCES = sim-main.c relrod.c utils.c image.c cell.c \
-                            hdf5-file.c
+                            hdf5-file.c templates.c
 simulate_patterns_LDADD = @LIBS@ -lm -lgsl -lgslcblas

src/hdf5-file.c

@@ -19,6 +19,8 @@
 #include <stdint.h>
 #include <hdf5.h>
 
+#include "image.h"
+
 
 int hdf5_write(const char *filename, const uint16_t *data,
                int width, int height)
@@ -74,3 +76,59 @@ int hdf5_write(const char *filename, const uint16_t *data,
 
 	return 0;
 }
+
+
+int hdf5_read(struct image *image, const char *filename)
+{
+	hid_t fh, sh, dh;	/* File, dataspace and data handles */
+	herr_t r;
+	hsize_t size[2];
+	hsize_t max_size[2];
+	uint16_t *buf;
+
+	fh = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+	if ( fh < 0 ) {
+		/* TODO: Try other formats here. */
+		fprintf(stderr, "Couldn't open file: %s\n", filename);
+		return 1;
+	}
+
+	dh = H5Dopen(fh, "/data/data", H5P_DEFAULT);
+	if ( dh < 0 ) {
+		fprintf(stderr, "Couldn't open dataset\n");
+		H5Fclose(fh);
+		return 1;
+	}
+
+	sh = H5Dget_space(dh);
+	if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
+		fprintf(stderr, "Dataset is not two-dimensional\n");
+		H5Fclose(fh);
+		return 1;
+	}
+
+	H5Sget_simple_extent_dims(sh, size, max_size);
+	printf("Data dimensions %i %i (max %i %i)\n",
+	                           (int)size[1], (int)size[0],
+	                           (int)max_size[1], (int)max_size[0]);
+
+	buf = malloc(sizeof(float)*size[0]*size[1]);
+
+	r = H5Dread(dh, H5T_NATIVE_UINT16, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf);
+	if ( r < 0 ) {
+		fprintf(stderr, "Couldn't read data\n");
+		H5Dclose(dh);
+		H5Fclose(fh);
+		return 1;
+	}
+
+	image->data = buf;
+	image->height = size[0];
+	image->width = size[1];
+	image->x_centre = image->width/2;
+	image->y_centre = image->height/2;
+
+	H5Fclose(fh);
+
+	return 0;
+}

src/hdf5-file.h

@@ -21,4 +21,6 @@
 extern int hdf5_write(const char *filename, const uint16_t *data,
                       int width, int height);
 
+extern int hdf5_read(struct image *image, const char *filename);
+
 #endif	/* HDF5_H */

src/main.c

@@ -20,6 +20,12 @@
 #include <string.h>
 #include <unistd.h>
 
+#include "cell.h"
+#include "image.h"
+#include "utils.h"
+#include "hdf5-file.h"
+#include "templates.h"
+
 
 static void main_show_help(const char *s)
 {
@@ -35,6 +41,9 @@ int main(int argc, char *argv[])
 	char **in_files;
 	size_t nin;
 	size_t i;
+	UnitCell *cell;
+	TemplateList *templates;
+	struct image template_parameters;
 
 	while ((c = getopt(argc, argv, "h")) != -1) {
 
@@ -60,11 +69,49 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
-	printf("Generating templates...\n");
+	/* Define unit cell */
+	cell = cell_new_from_parameters(28.10e-9,
+	                                28.10e-9,
+	                                16.52e-9,
+	                                deg2rad(90.0),
+	                                deg2rad(90.0),
+	                                deg2rad(120.0));
+
+	/* Generate templates */
+	template_parameters.width = 512;
+	template_parameters.height = 512;
+	template_parameters.fmode = FORMULATION_CLEN;
+	template_parameters.x_centre = 255.5;
+	template_parameters.y_centre = 255.5;
+	template_parameters.camera_len = 0.2;  /* 20 cm */
+	template_parameters.resolution = 5120; /* 512 pixels in 10 cm */
+	template_parameters.lambda = 0.2e-9;   /* LCLS wavelength */
+	templates = generate_templates(cell, template_parameters);
 
 	printf("Input files (%i):\n", nin);
 	for ( i=0; i<nin; i++ ) {
+
+		struct image image;
+
 		printf("%6i: %s\n", i+1, in_files[i]);
+
+		image.width = 512;
+		image.height = 512;
+		image.fmode = FORMULATION_CLEN;
+		image.x_centre = 255.5;
+		image.y_centre = 255.5;
+		image.camera_len = 0.2;  /* 20 cm */
+		image.resolution = 5120; /* 512 pixels in 10 cm */
+		image.lambda = 0.2e-9;   /* LCLS wavelength */
+
+		if ( hdf5_read(&image, in_files[i]) ) {
+			fprintf(stderr, "Couldn't read file '%s'\n",
+			        in_files[i]);
+			continue;
+		}
+
+		try_templates(&image, templates);
+
 	}
 
 	return 0;

src/templates.c

+/*
+ * templates.c
+ *
+ * Handle templates
+ *
+ * (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 <stdio.h>
+
+#include "templates.h"
+#include "cell.h"
+#include "image.h"
+#include "utils.h"
+#include "relrod.h"
+
+
+struct _templatelist
+{
+	int		n_templates;
+	struct template	*templates;
+};
+
+
+struct template_feature
+{
+	float		x;
+	float		y;
+	struct template_feature *next;
+};
+
+
+struct template
+{
+	float		omega;
+	float		tilt;
+	struct template_feature *features;
+};
+
+
+static int template_add(TemplateList *list, struct template *template)
+{
+	if ( list->templates ) {
+		list->templates = realloc(list->templates,
+		                 (list->n_templates+1)*sizeof(struct template));
+	} else {
+		assert(list->n_templates == 0);
+		list->templates = malloc(sizeof(struct template));
+	}
+
+	/* Copy the data */
+	list->templates[list->n_templates] = *template;
+
+	list->n_templates++;
+
+	return list->n_templates - 1;
+}
+
+
+static TemplateList *template_list_new()
+{
+	TemplateList *list;
+
+	list = malloc(sizeof(TemplateList));
+
+	list->n_templates = 0;
+	list->templates = NULL;
+
+	return list;
+}
+
+
+TemplateList *generate_templates(UnitCell *cell, struct image params)
+{
+	TemplateList *list;
+	double omega, tilt;
+
+	list = template_list_new();
+
+	omega = deg2rad(40.0);
+
+	//for ( omega=deg2rad(-180); omega<deg2rad(180); omega+=deg2rad(1) ) {
+
+		params.omega = omega;
+
+		for ( tilt=0; tilt<deg2rad(180); tilt+=deg2rad(1) ) {
+
+			struct template t;
+			struct template_feature *tfc;
+			int nrefl, i;
+
+			t.omega = omega;
+			t.tilt = tilt;
+			t.features = malloc(sizeof(struct template_feature));
+			t.features->next = NULL;
+			tfc = t.features;
+
+			params.tilt = tilt;
+			get_reflections(&params, cell, 0.01e9);
+
+			nrefl = image_feature_count(params.rflist);
+			for ( i=0; i<nrefl; i++ ) {
+
+				struct imagefeature *f;
+
+				f = image_get_feature(params.rflist, i);
+
+				tfc->x = f->x;
+				tfc->y = f->y;
+				tfc->next = malloc(sizeof(struct template_feature));
+				tfc = tfc->next;
+
+			}
+
+			template_add(list, &t);
+
+			image_feature_list_free(params.rflist);
+
+			printf("Generating templates... %+5.2f %+5.2f\r",
+			       rad2deg(omega), rad2deg(tilt));
+
+		}
+	//}
+
+	return list;
+}
+
+
+int try_templates(struct image *image, TemplateList *list)
+{
+	return 0;
+}

src/templates.h

+/*
+ * templates.h
+ *
+ * Handle templates
+ *
+ * (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 TEMPLATES_H
+#define TEMPLATES_H
+
+#include "cell.h"
+#include "image.h"
+
+/* An opaque type representing a list of templates */
+typedef struct _templatelist TemplateList;
+
+extern TemplateList *generate_templates(UnitCell *cell, struct image params);
+extern int try_templates(struct image *image, TemplateList *list);
+
+#endif	/* TEMPLAETS_H */