Add intensity output

4 files changed, 168 insertions, 1 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 7a8e9e43..47f6f002 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -3,7 +3,7 @@ bin_PROGRAMS = pattern_sim integr_sim
 AM_CFLAGS = -Wall
 
 pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c cell.c \
-                      hdf5-file.c ewald.c detector.c sfac.c
+                      hdf5-file.c ewald.c detector.c sfac.c intensities.c
 pattern_sim_LDADD = @LIBS@
 
 integr_sim_SOURCES = integr_sim.c diffraction.c utils.c ewald.c sfac.c \
diff --git a/src/intensities.c b/src/intensities.c
new file mode 100644
index 00000000..7630bcd8
--- /dev/null
+++ b/src/intensities.c
@@ -0,0 +1,139 @@
+/*
+ * intensities.c
+ *
+ * Extract intensities from patterns
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * pattern_sim - Simulate diffraction patterns from small crystals
+ *
+ */
+
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <assert.h>
+
+#include "image.h"
+#include "intensities.h"
+#include "cell.h"
+#include "sfac.h"
+
+
+#define MAX_HITS (1024)
+
+
+struct reflhit {
+	signed int h;
+	signed int k;
+	signed int l;
+	double min_distance;
+	int x;
+	int y;
+};
+
+
+static int sum_nearby_points(uint16_t *data, int width, int x, int y)
+{
+	int dx, dy;
+	int intensity = 0;
+
+	for ( dx=-3; dx<=3; dx++ ) {
+	for ( dy=-3; dy<=3; dy++ ) {
+		intensity += data[(x+dx) + width*(y+dy)];
+		data[(x+dx) + width*(y+dy)] = 0;
+	}
+	}
+
+	return intensity;
+}
+
+
+void output_intensities(struct image *image)
+{
+	int x, y;
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+	struct reflhit hits[MAX_HITS];
+	int n_hits = 0;
+	int i;
+
+	cell_get_cartesian(image->molecule->cell, &ax, &ay, &az,
+		                                  &bx, &by, &bz,
+		                                  &cx, &cy, &cz);
+
+	for ( x=0; x<image->width; x++ ) {
+	for ( y=0; y<image->height; y++ ) {
+
+		double hd, kd, ld;
+		signed int h, k, l;
+		struct threevec q;
+		double dist;
+		int found = 0;
+		int j;
+
+		q = image->qvecs[x + image->width*y];
+
+		hd = q.u * ax + q.v * ay + q.w * az;
+		kd = q.u * bx + q.v * by + q.w * bz;
+		ld = q.u * cx + q.v * cy + q.w * cz;
+
+		h = (signed int)rint(hd);
+		k = (signed int)rint(kd);
+		l = (signed int)rint(ld);
+
+		dist = sqrt(pow(fmod(hd, 1.0), 2.0) + pow(fmod(kd, 1.0), 2.0)
+		          + pow(fmod(ld, 1.0), 2.0));
+
+		if ( dist > 0.2 ) continue;
+
+		for ( j=0; j<n_hits; j++ ) {
+			if ( (hits[j].h == h) && (hits[j].k == k)
+			                      && (hits[j].l == l) ) {
+
+				if ( dist < hits[j].min_distance ) {
+					hits[j].min_distance = dist;
+					hits[j].x = x;
+					hits[j].y = y;
+				}
+
+				found = 1;
+
+			}
+		}
+
+		if ( !found ) {
+			hits[n_hits].min_distance = dist;
+			hits[n_hits].x = x;
+			hits[n_hits].y = y;
+			hits[n_hits].h = h;
+			hits[n_hits].k = k;
+			hits[n_hits].l = l;
+			n_hits++;
+			assert(n_hits < MAX_HITS);
+		}
+
+	}
+	}
+
+	STATUS("Found %i reflections\n", n_hits);
+
+	/* Explicit printf() used here (not normally allowed) because
+	 * we really want to output to stdout */
+	for ( i=0; i<n_hits; i++ ) {
+
+		int intensity;
+		intensity = sum_nearby_points(image->data, image->width,
+		                              hits[i].x, hits[i].y);
+
+		printf("%3i %3i %3i %6i (at %i,%i)\n",
+		       hits[i].h, hits[i].k, hits[i].l, intensity,
+		       hits[i].x, hits[i].y);
+
+	}
+
+	/* Blank line at end */
+	printf("\n");
+}
diff --git a/src/intensities.h b/src/intensities.h
new file mode 100644
index 00000000..1a788f53
--- /dev/null
+++ b/src/intensities.h
@@ -0,0 +1,23 @@
+/*
+ * intensities.h
+ *
+ * Extract intensities from patterns
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * pattern_sim - Simulate diffraction patterns from small crystals
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef INTENSITIES_H
+#define INTENSITIES_H
+
+#include "image.h"
+
+extern void output_intensities(struct image *image);
+
+#endif	/* INTENSITIES_H */
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 3f3d93ea..e3c5ec8a 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -25,6 +25,7 @@
 #include "utils.h"
 #include "hdf5-file.h"
 #include "detector.h"
+#include "intensities.h"
 
 
 static void show_help(const char *s)
@@ -224,6 +225,10 @@ int main(int argc, char *argv[])
 		get_diffraction(&image);
 		record_image(&image);
 
+		if ( config_nearbragg ) {
+			output_intensities(&image);
+		}
+
 		snprintf(filename, 1023, "results/sim-%i.h5", number);
 		number++;