Rework multisampling (remove Ewald module)

10 files changed, 81 insertions, 267 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index a80ba7a2..6a40c9a0 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -8,7 +8,7 @@ AM_CFLAGS = -Wall
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\"
 
 pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c cell.c \
-                      hdf5-file.c ewald.c detector.c sfac.c intensities.c \
+                      hdf5-file.c detector.c sfac.c intensities.c \
                       reflections.c
 pattern_sim_LDADD = @LIBS@
 
@@ -17,7 +17,7 @@ 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 dirax.c cell.c image.c \
-                      intensities.c ewald.c peaks.c index.c filters.c \
+                      intensities.c peaks.c index.c filters.c \
                       diffraction.c detector.c sfac.c
 indexamajig_LDADD = @LIBS@
 
diff --git a/src/detector.c b/src/detector.c
index ccf37024..ae7576b2 100644
--- a/src/detector.c
+++ b/src/detector.c
@@ -182,7 +182,7 @@ void record_image(struct image *image, int do_water, int do_poisson,
 		if ( do_water ) {
 
 			/* Add intensity contribution from water */
-			water = water_intensity(image->qvecs[0][x + image->width*y],
+			water = water_intensity(get_q(image, x, y, 1, NULL),
 			                        ph_lambda_to_en(image->lambda),
 			                        BEAM_RADIUS, WATER_RADIUS);
 			intensity += water;
diff --git a/src/diffraction.c b/src/diffraction.c
index b3d6bef0..bc3e685c 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -19,15 +19,17 @@
 #include "image.h"
 #include "utils.h"
 #include "cell.h"
-#include "ewald.h"
 #include "diffraction.h"
 #include "sfac.h"
 
 
+#define SAMPLING (1)
+
+
 static double lattice_factor(struct rvec q, double ax, double ay, double az,
-                                                double bx, double by, double bz,
-                                                double cx, double cy, double cz,
-                                                int na, int nb, int nc)
+                                            double bx, double by, double bz,
+                                            double cx, double cy, double cz,
+                                            int na, int nb, int nc)
 {
 	struct rvec Udotq;
 	double f1, f2, f3;
@@ -133,19 +135,56 @@ double water_intensity(struct rvec q, double en,
 }
 
 
+struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
+                  unsigned int sampling, float *ttp)
+{
+	struct rvec q;
+	float twothetax, twothetay, twotheta, r;
+	float rx = 0.0;
+	float ry = 0.0;
+	int p;
+
+	const float k = 1.0/image->lambda;
+	const unsigned int x = xs / sampling;
+	const unsigned int y = ys / sampling; /* Integer part only */
+
+	for ( p=0; p<image->det.n_panels; p++ ) {
+		if ( (x >= image->det.panels[p].min_x)
+		  && (x <= image->det.panels[p].max_x)
+		  && (y >= image->det.panels[p].min_y)
+		  && (y <= image->det.panels[p].max_y) ) {
+			rx = ((float)xs - (sampling*image->det.panels[p].cx))
+			               / (sampling * image->resolution);
+			ry = ((float)ys - (sampling*image->det.panels[p].cy))
+			               / (sampling * image->resolution);
+			break;
+		}
+	}
+
+	/* Calculate q-vector for this sub-pixel */
+	r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
+	twothetax = atan2(rx, image->camera_len);
+	twothetay = atan2(ry, image->camera_len);
+	twotheta = atan2(r, image->camera_len);
+
+	if ( ttp != NULL ) *ttp = twotheta;
+
+	q.u = k * sin(twothetax);
+	q.v = k * sin(twothetay);
+	q.w = k - k * cos(twotheta);
+
+	return q;
+}
+
+
 void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac)
 {
-	int x, y;
+	unsigned int xs, ys;
 	double ax, ay, az;
 	double bx, by, bz;
 	double cx, cy, cz;
 	double a, b, c, d;
 
-	/* Generate the array of reciprocal space vectors in image->qvecs */
-	if ( image->qvecs == NULL ) {
-		get_ewald(image);
-	}
-
 	if ( image->molecule == NULL ) return;
 
 	cell_get_cartesian(image->molecule->cell, &ax, &ay, &az,
@@ -157,6 +196,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac)
 	STATUS("Particle size = %i x %i x %i (=%5.2f x %5.2f x %5.2f nm)\n",
 	       na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9);
 
+	/* Allocate (and zero) the "diffraction array" */
 	image->sfacs = calloc(image->width * image->height,
 	                      sizeof(double complex));
 
@@ -167,34 +207,36 @@ void get_diffraction(struct image *image, int na, int nb, int nc, int no_sfac)
 		}
 	}
 
-	for ( x=0; x<image->width; x++ ) {
-	for ( y=0; y<image->height; y++ ) {
+	/* Needed later for Lorentz calculation */
+	image->twotheta = malloc(image->width * image->height * sizeof(double));
+
+	for ( xs=0; xs<image->width*SAMPLING; xs++ ) {
+	for ( ys=0; ys<image->height*SAMPLING; ys++ ) {
 
 		double f_lattice;
 		double complex f_molecule;
 		struct rvec q;
-		double complex val;
-		int i;
-
-		for ( i=0; i<image->nspheres; i++ ) {
-
-			q = image->qvecs[i][x + image->width*y];
-
-			f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz,
-			                           na, nb, nc);
-			if ( no_sfac ) {
-				f_molecule = 10000.0;
-			} else {
-				f_molecule = molecule_factor(image->molecule, q,
-			                            ax,ay,az,bx,by,bz,cx,cy,cz);
-			}
-
-			val = f_molecule * f_lattice;
-			image->sfacs[x + image->width*y] += val;
-
+		float twotheta;
+		double sw = 1.0/(SAMPLING*SAMPLING);  /* Sample weight */
+
+		const unsigned int x = xs / SAMPLING;
+		const unsigned int y = ys / SAMPLING; /* Integer part only */
+
+		q = get_q(image, xs, ys, SAMPLING, &twotheta);
+		image->twotheta[x + image->width*y] = twotheta;
+
+		f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz,
+		                           na, nb, nc);
+		if ( no_sfac ) {
+			f_molecule = 10000.0;
+		} else {
+			f_molecule = molecule_factor(image->molecule, q,
+		                            ax,ay,az,bx,by,bz,cx,cy,cz);
 		}
 
+		image->sfacs[x + image->width*y] += sw * f_molecule * f_lattice;
+
 	}
-	progress_bar(x, image->width-1, "Calculating lattice factors");
+	progress_bar(xs, SAMPLING*image->width-1, "Calculating lattice factors");
 	}
 }
diff --git a/src/diffraction.h b/src/diffraction.h
index efc59d36..36d1c21f 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -23,5 +23,6 @@ extern void get_diffraction(struct image *image, int na, int nb, int nc,
                             int nosfac);
 extern double water_intensity(struct rvec q, double en,
                               double beam_r, double water_r);
-
+extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
+                         unsigned int sampling, float *ttp);
 #endif	/* DIFFRACTION_H */
diff --git a/src/ewald.c b/src/ewald.c
deleted file mode 100644
index fca488f6..00000000
--- a/src/ewald.c
+++ /dev/null
@@ -1,196 +0,0 @@
-/*
- * ewald.c
- *
- * Calculate q-vector arrays
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-
-#include "image.h"
-#include "utils.h"
-#include "cell.h"
-#include "ewald.h"
-#include "detector.h"
-
-
-#define SAMPLING (4)
-#define BWSAMPLING (10)
-#define BANDWIDTH (0.05)
-
-
-static struct rvec quat_rot(struct rvec q, struct quaternion z)
-{
-	struct rvec res;
-	double t01, t02, t03, t11, t12, t13, t22, t23, t33;
-
-	t01 = z.w*z.x;
-	t02 = z.w*z.y;
-	t03 = z.w*z.z;
-	t11 = z.x*z.x;
-	t12 = z.x*z.y;
-	t13 = z.x*z.z;
-	t22 = z.y*z.y;
-	t23 = z.y*z.z;
-	t33 = z.z*z.z;
-
-	res.u = (1.0 - 2.0 * (t22 + t33)) * q.u
-	            + (2.0 * (t12 + t03)) * q.v
-	            + (2.0 * (t13 - t02)) * q.w;
-
-	res.v =       (2.0 * (t12 - t03)) * q.u
-	      + (1.0 - 2.0 * (t11 + t33)) * q.v
-	            + (2.0 * (t01 + t23)) * q.w;
-
-	res.w =       (2.0 * (t02 + t13)) * q.u
-	            + (2.0 * (t23 - t01)) * q.v
-	      + (1.0 - 2.0 * (t11 + t22)) * q.w;
-
-	return res;
-}
-
-
-static void add_sphere(struct image *image, double k, int soffs)
-{
-	int x, y;
-
-	for ( x=0; x<image->width; x++ ) {
-	for ( y=0; y<image->height; y++ ) {
-
-		double rx = 0.0;
-		double ry = 0.0;
-		double r;
-		double twothetax, twothetay, twotheta;
-		double qx, qy, qz;
-		struct rvec q1, q2, q3, q4;
-		int p, sx, sy, i;
-
-		/* Calculate q vectors for Ewald sphere */
-		for ( p=0; p<image->det.n_panels; p++ ) {
-			if ( (x >= image->det.panels[p].min_x)
-			  && (x <= image->det.panels[p].max_x)
-			  && (y >= image->det.panels[p].min_y)
-			  && (y <= image->det.panels[p].max_y) ) {
-				rx = ((double)x - image->det.panels[p].cx)
-				                            / image->resolution;
-				ry = ((double)y - image->det.panels[p].cy)
-				                            / image->resolution;
-				break;
-			}
-		}
-
-		/* Bottom left corner */
-		r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
-		twothetax = atan2(rx, image->camera_len);
-		twothetay = atan2(ry, image->camera_len);
-		twotheta = atan2(r, image->camera_len);
-		qx = k * sin(twothetax);
-		qy = k * sin(twothetay);
-		qz = k - k * cos(twotheta);
-		q1.u = qx;  q1.v = qy;  q1.w = qz;
-		/* 2theta value is calculated at the bottom left only */
-		image->twotheta[x + image->width*y] = twotheta;
-
-		/* Bottom right corner (using the same panel configuration!) */
-		rx = ((double)(x+1) - image->det.panels[p].cx)
-		                            / image->resolution;
-		ry = ((double)y - image->det.panels[p].cy)
-		                            / image->resolution;
-		twothetax = atan2(rx, image->camera_len);
-		twothetay = atan2(ry, image->camera_len);
-		twotheta = atan2(r, image->camera_len);
-		qx = k * sin(twothetax);
-		qy = k * sin(twothetay);
-		qz = k - k * cos(twotheta);
-		q2.u = qx;  q2.v = qy;  q2.w = qz;
-
-		/* Top left corner (using the same panel configuration!) */
-		rx = ((double)x - image->det.panels[p].cx)
-		                            / image->resolution;
-		ry = ((double)(y+1) - image->det.panels[p].cy)
-		                            / image->resolution;
-		twothetax = atan2(rx, image->camera_len);
-		twothetay = atan2(ry, image->camera_len);
-		twotheta = atan2(r, image->camera_len);
-		qx = k * sin(twothetax);
-		qy = k * sin(twothetay);
-		qz = k - k * cos(twotheta);
-		q3.u = qx;  q3.v = qy;  q3.w = qz;
-
-		/* Top right corner (using the same panel configuration!) */
-		rx = ((double)(x+1) - image->det.panels[p].cx)
-		                            / image->resolution;
-		ry = ((double)(y+1) - image->det.panels[p].cy)
-		                            / image->resolution;
-		twothetax = atan2(rx, image->camera_len);
-		twothetay = atan2(ry, image->camera_len);
-		twotheta = atan2(r, image->camera_len);
-		qx = k * sin(twothetax);
-		qy = k * sin(twothetay);
-		qz = k - k * cos(twotheta);
-		q4.u = qx;  q4.v = qy;  q4.w = qz;
-
-		/* Now interpolate between the values to get
-		 * the sampling points */
-		i = soffs;
-		for ( sx=0; sx<SAMPLING; sx++ ) {
-		for ( sy=0; sy<SAMPLING; sy++ ) {
-
-			struct rvec q;
-
-			q.u = q1.u + ((q2.u - q1.u)/SAMPLING)*sx
-			           + ((q3.u - q1.u)/SAMPLING)*sy;
-			q.v = q1.v + ((q2.v - q1.v)/SAMPLING)*sx
-			           + ((q3.v - q1.v)/SAMPLING)*sy;
-			q.w = q1.w + ((q2.w - q1.w)/SAMPLING)*sx
-			           + ((q3.w - q1.w)/SAMPLING)*sy;
-			image->qvecs[i++][x + image->width*y] = quat_rot(q,
-		                                            image->orientation);
-
-		}
-		}
-
-	}
-	}
-
-}
-
-
-void get_ewald(struct image *image)
-{
-	double kc;  /* Wavenumber */
-	int i, kstep;
-	long long int mtotal = 0;
-
-	kc = 1/image->lambda;  /* Centre */
-
-	image->twotheta = malloc(image->width * image->height
-	                       * sizeof(double));
-
-	/* Create the spheres */
-	image->nspheres = SAMPLING*SAMPLING*BWSAMPLING;
-	image->qvecs = malloc(image->nspheres * sizeof(struct rvec *));
-
-	for ( i=0; i<image->nspheres; i++ ) {
-		mtotal += image->width * image->height * sizeof(struct rvec);
-		image->qvecs[i] = malloc(image->width * image->height
-                                                 * sizeof(struct rvec));
-	}
-	STATUS("%i spheres, %lli Mbytes\n", image->nspheres, mtotal/(1024*1024));
-
-	for ( kstep=0; kstep<BWSAMPLING; kstep++ ) {
-
-		double k;
-
-		k = kc + (kstep-(BWSAMPLING/2))*kc*(BANDWIDTH/BWSAMPLING);
-		add_sphere(image, k, kstep*SAMPLING);
-
-	}
-}
diff --git a/src/ewald.h b/src/ewald.h
deleted file mode 100644
index 75ebb44d..00000000
--- a/src/ewald.h
+++ /dev/null
@@ -1,23 +0,0 @@
-/*
- * ewald.h
- *
- * Calculate q-vector arrays
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef EWALD_H
-#define EWALD_H
-
-#include "image.h"
-
-extern void get_ewald(struct image *image);
-
-#endif	/* EWALD_H */
diff --git a/src/image.h b/src/image.h
index ff87b1bc..144ebe42 100644
--- a/src/image.h
+++ b/src/image.h
@@ -76,8 +76,6 @@ struct image {
 	int			*hdr;      /* Actual counts */
 	int16_t			*data;     /* Integer counts after bloom */
 	double complex		*sfacs;
-	struct rvec		**qvecs;
-	int			nspheres;
 	double			*twotheta;
 	struct molecule		*molecule;
 	UnitCell		*indexed_cell;
diff --git a/src/indexamajig.c b/src/indexamajig.c
index d05db543..43127e66 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -26,7 +26,6 @@
 #include "hdf5-file.h"
 #include "index.h"
 #include "intensities.h"
-#include "ewald.h"
 #include "peaks.h"
 #include "diffraction.h"
 #include "detector.h"
@@ -236,7 +235,6 @@ int main(int argc, char *argv[])
 
 				/* Simulate a diffraction pattern */
 				image.sfacs = NULL;
-				image.qvecs = NULL;
 				image.twotheta = NULL;
 				image.hdr = NULL;
 
@@ -245,7 +243,6 @@ int main(int argc, char *argv[])
 				image.orientation.x = 0.0;
 				image.orientation.y = 0.0;
 				image.orientation.z = 0.0;
-				get_ewald(&image);
 
 			}
 
diff --git a/src/intensities.c b/src/intensities.c
index edc15910..8b5f51c5 100644
--- a/src/intensities.c
+++ b/src/intensities.c
@@ -19,6 +19,7 @@
 #include "intensities.h"
 #include "cell.h"
 #include "sfac.h"
+#include "diffraction.h"
 
 
 #define MAX_HITS (1024)
@@ -72,7 +73,7 @@ void output_intensities(struct image *image, UnitCell *cell)
 		int found = 0;
 		int j;
 
-		q = image->qvecs[0][x + image->width*y];
+		q = get_q(image, x, y, 1, NULL);
 
 		hd = q.u * ax + q.v * ay + q.w * az;
 		kd = q.u * bx + q.v * by + q.w * bz;
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 70384250..5fd71fc1 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -247,7 +247,6 @@ int main(int argc, char *argv[])
 	do {
 
 		int na, nb, nc;
-		int i;
 
 		na = 8*random()/RAND_MAX + 4;
 		nb = 8*random()/RAND_MAX + 4;
@@ -270,7 +269,6 @@ int main(int argc, char *argv[])
 		}
 
 		/* Ensure no residual information */
-		image.qvecs = NULL;
 		image.sfacs = NULL;
 		image.data = NULL;
 		image.twotheta = NULL;
@@ -303,10 +301,6 @@ int main(int argc, char *argv[])
 
 		/* Clean up */
 		free(image.data);
-		for ( i=0; i<image.nspheres; i++ ) {
-			free(image.qvecs[i]);
-		}
-		free(image.qvecs);
 		free(image.hdr);
 		free(image.sfacs);
 		free(image.twotheta);