Implement GRADIENT_MOSAIC and GRADIENT_INTERPOLATE

1 files changed, 103 insertions, 11 deletions

diff --git a/src/diffraction.c b/src/diffraction.c
index dcd6ae2c..1960f786 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -15,6 +15,7 @@
 #include <stdio.h>
 #include <string.h>
 #include <complex.h>
+#include <assert.h>
 
 #include "image.h"
 #include "utils.h"
@@ -69,31 +70,121 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az,
 }
 
 
+static double interpolate_linear(const double *ref,
+                                 const unsigned int *counts,
+                                 float hd, signed int k, signed int l)
+{
+	signed int h;
+	double val1, val2;
+	float f;
+	unsigned int c1, c2;
+
+	h = (signed int)hd;
+	if ( hd < 0.0 ) h -= 1;
+	f = hd - (float)h;
+	assert(f >= 0.0);
+
+	val1 = lookup_intensity(ref, h, k, l);
+	val2 = lookup_intensity(ref, h+1, k, l);
+	c1 = lookup_count(counts, h, k, l);
+	c2 = lookup_count(counts, h+1, k, l);
+
+	if ( c1 == 0 ) {
+		ERROR("Needed intensity for %i %i %i, but don't have it.\n",
+		      h, k, l);
+		return 1.0e20;
+	}
+
+	if ( c2 == 0 ) {
+		ERROR("Needed intensity for %i %i %i, but don't have it.\n",
+		      h+1, k, l);
+		return 1.0e20;
+	}
+
+	val1 = val1 / (double)c1;
+	val2 = val2 / (double)c2;
+
+	return (1.0-f)*val1 + f*val2;
+}
+
+
+static double interpolate_bilinear(const double *ref,
+                                   const unsigned int *counts,
+                                   float hd, float kd, signed int l)
+{
+	signed int k;
+	double val1, val2;
+	float f;
+
+	k = (signed int)kd;
+	if ( kd < 0.0 ) k -= 1;
+	f = kd - (float)k;
+	assert(f >= 0.0);
+
+	val1 = interpolate_linear(ref, counts, hd, k, l);
+	val2 = interpolate_linear(ref, counts, hd, k+1, l);
+
+	return (1.0-f)*val1 + f*val2;
+}
+
+
+static double interpolate_intensity(const double *ref,
+                                    const unsigned int *counts,
+                                    float hd, float kd, float ld)
+{
+	signed int l;
+	double val1, val2;
+	float f;
+
+	l = (signed int)ld;
+	if ( ld < 0.0 ) l -= 1;
+	f = ld - (float)l;
+	assert(f >= 0.0);
+
+	val1 = interpolate_bilinear(ref, counts, hd, kd, l);
+	val2 = interpolate_bilinear(ref, counts, hd, kd, l+1);
+
+	return (1.0-f)*val1 + f*val2;
+}
+
+
 /* Look up the structure factor for the nearest Bragg condition */
 static double molecule_factor(const double *intensities,
                               const unsigned int *counts, struct rvec q,
                               double ax, double ay, double az,
                               double bx, double by, double bz,
-                              double cx, double cy, double cz)
+                              double cx, double cy, double cz,
+                              GradientMethod m)
 {
-	double hd, kd, ld;
+	float hd, kd, ld;
 	signed int h, k, l;
 	double r;
 
 	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);
 
-	if ( lookup_count(counts, h, k, l) == 0 ) {
-		ERROR("Needed intensity for %i %i %i, but don't have it.\n",
-		      h, k, l);
-		return 1.0e20;
+	switch ( m ) {
+	case GRADIENT_MOSAIC :
+		h = (signed int)rint(hd);
+		k = (signed int)rint(kd);
+		l = (signed int)rint(ld);
+		if ( lookup_count(counts, h, k, l) == 0 ) {
+			ERROR("Needed intensity for %i %i %i, but don't have it.\n",
+			      h, k, l);
+			return 1.0e20;
+		}
+		r = lookup_intensity(intensities, h, k, l);
+		break;
+	case GRADIENT_INTERPOLATE :
+		r = interpolate_intensity(intensities, counts, hd, kd, ld);
+		break;
+	case GRADIENT_PHASED :
+	default:
+		ERROR("This gradient method not implemented yet.\n");
+		exit(1);
 	}
 
-	r = lookup_intensity(intensities, h, k, l);
 	return r;
 }
 
@@ -235,7 +326,8 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
 				I_molecule = molecule_factor(intensities,
 				                             counts, q,
 			                                     ax,ay,az,
-			                                     bx,by,bz,cx,cy,cz);
+			                                     bx,by,bz,cx,cy,cz,
+				                             m);
 			}
 
 			I_lattice = pow(f_lattice, 2.0);