More de-clobbering of reprojection

1 files changed, 132 insertions, 74 deletions

diff --git a/src/reproject.c b/src/reproject.c
index df85f4f..e0199af 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -18,6 +18,7 @@
 #include "imagedisplay.h"
 #include "displaywindow.h"
 #include "image.h"
+#include "mapping.h"
 
 /* Attempt to find partners from the feature list of "image" for each feature in "flist". */
 void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) {
@@ -37,6 +38,8 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) {
 		if ( (d <= 20.0) && partner ) {
 			rflist->features[i].partner = partner;
 			rflist->features[i].partner_d = d;
+			image->features->features[idx].partner =
+				&rflist->features[i];
 		} else {
 			rflist->features[i].partner = NULL;
 		}
@@ -45,16 +48,19 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) {
 
 }
 
-ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist) {
-
+ImageFeatureList *reproject_get_reflections(ImageRecord *image,
+						ReflectionList *reflectionlist)
+{
 	ImageFeatureList *flist;
 	Reflection *reflection;
-	double smax = 0.4e9;
+	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, uxt, uyt, uzt;	/* "up" vector (and calculation intermediates) */
-	//int done_debug = 0;
+	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 */
+	int nrr = 0;
+	int n = 0;
 
 	flist = image_feature_list_new();
 
@@ -67,15 +73,12 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 	kx = nx / image->lambda;
 	ky = ny / image->lambda;
 	kz = nz / image->lambda;	/* This is the centre of the Ewald sphere */
-	//reflection_add(ctx->reflectionlist, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1);
 
-	/* Determine where "up" is
-	 *	See above. */
-	uxt = 0.0;  uyt = 1.0;  uzt = 0.0;
-	ux = uxt;  uy = cos(tilt)*uyt + sin(tilt)*uzt;  uz = -sin(tilt)*uyt + cos(tilt)*uzt;
-	uxt = ux;  uyt = uy;  uzt = uz;
-	ux = uxt*cos(-omega) + uyt*-sin(omega);  uy = -uxt*sin(omega) + uyt*cos(omega);  uz = uzt;
-	//reflection_add(ctx->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2);
+	/* 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);
 
 	reflection = reflectionlist->reflections;
 	do {
@@ -92,7 +95,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 		g_sq = modulus_squared(xl, yl, zl);
 		gn = xl*nx + yl*ny + zl*nz;
 
-		/* Next, solve the relrod equation to calculate the excitation error */
+		/* 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;
@@ -106,10 +110,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 
 			double xi, yi, zi;
 			double gx, gy, gz;
-			double cx, cy, cz;
 			double theta;
 			double x, y;
-			double rx, ry, rz;
 
 			/* Determine the intersection point */
 			xi = xl + s*nx;  yi = yl + s*ny;  zi = zl + s*nz;
@@ -117,64 +119,36 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 			/* 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 */
+			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 psi, disc;
-
-				//reflection_add(ctx->reflectionlist, xl, yl, zl, 1, REFLECTION_GENERATED);
-				//reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_MARKER);
-
-				/* Calculate azimuth of point in image (clockwise from "up", will be changed later) */
-				cx = yi*nz-zi*ny;  cy = nx*zi-nz*xi;  cz = ny*xi-nx*yi;  /* c = i x n */
-				psi = angle_between(cx, cy, cz, ux, uy, uz);
-
-				/* Finally, resolve the +/- Pi ambiguity from the previous step */
-				rx = cy*nz-cz*ny;  ry = nx*cz-nz*cx;  rz = ny*cx-nx*cy;  /* r = [i x n] x n */
-				disc = angle_between(rx, ry, rz, ux, uy, uz);
-			//	if ( (i==20) && !done_debug ) {
-			//		reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_VECTOR_MARKER_3);
-			//		reflection_add(ctx->reflectionlist, cx, cy, cz, 1, REFLECTION_VECTOR_MARKER_4);
-			//		reflection_add(ctx->reflectionlist, rx, ry, rz, 1, REFLECTION_VECTOR_MARKER_4);
-			//		printf("psi=%f deg, disc=%f deg\n", rad2deg(psi), rad2deg(disc));
-			//	}
-				if ( (psi >= M_PI_2) && (disc >= M_PI_2) ) {
-					psi -= M_PI_2;		/* Case 1 */
-			//		if ( (i==20) && !done_debug ) printf("case 1\n");
-				} else if ( (psi >= M_PI_2) && (disc < M_PI_2) ) {
-					psi = 3*M_PI_2 - psi;	/* Case 2 */
-			//		if ( (i==20) && !done_debug ) printf("case 2\n");
-				} else if ( (psi < M_PI_2) && (disc < M_PI_2) ) {
-					psi = 3*M_PI_2 - psi;	/* Case 3 */
-			//		if ( (i==20) && !done_debug ) printf("case 3\n");
-				} else if ( (psi < M_PI_2) && (disc >= M_PI_2) ) {
-					psi = 3*M_PI_2 + psi;	/* Case 4 */
-			//		if ( (i==20) && !done_debug ) printf("case 4\n");
-				}
+				double dx, dy, psi;
 
-			//	if ( (i==20) && !done_debug ) printf("final psi=%f clockwise from 'up'\n", rad2deg(psi));
-				psi = M_PI_2 - psi;	/* Anticlockwise from "+x" instead of clockwise from "up" */
-			//	if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +x\n", rad2deg(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);
 
-				psi += omega;
-			//	if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +tilt axis\n", rad2deg(psi));
-			//	if ( (i==20) && !done_debug ) done_debug = 1;
-
-				/* Calculate image coordinates from polar representation */
+				/* 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 ) {
+				} 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");
+					fprintf(stderr,
+						"Unrecognised formulation mode "
+						"in reproject_get_reflections\n");
 					return NULL;
 				}
 
@@ -182,31 +156,45 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 				y += image->y_centre;
 
 				/* Sanity check */
-				if ( (x>=0) && (x<image->width) && (y>=0) && (y<image->height) ) {
-
-					/* Record the reflection */
-					image_add_feature_reflection(flist, x, y, image, reflection->intensity, reflection);
-					/* Intensity should be multiplied by relrod spike function except
-					 * reprojected reflections aren't used quantitatively for anything */
-
-					//printf("Reflection at %f, %f\n", x, y);
+				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, reflection->intensity,
+						reflection);
+					printf("adding %i %i %i\n",
+						reflection->h,
+						reflection->k,
+						reflection->l);
+					n++;
 
 				} /* else it's outside the picture somewhere */
 
-			} /* else this is the central beam so don't worry about it */
+			} /* else this is the central beam */
 
 		}
 
 		reflection = reflection->next;
 
+		nrr++;
+
 	} while ( reflection );
 
-	/* Partner features only if the image has a feature list.  This allows the test
-	 *	program to use this function to generate simulated data. */
+	/* Partner features only if the image has a feature list.  This allows
+	 *  the test programs to use this function to generate simulated data
+	 */
 	if ( image->features != NULL ) {
 		reproject_partner_features(flist, image);
 	}
 
+	printf("processed %i, found %i\n", nrr, n);
+
 	return flist;
 
 }
@@ -226,8 +214,20 @@ void reproject_cell_to_lattice(ControlContext *ctx) {
 
 	/* Invalidate all the reprojected feature lists */
 	for ( i=0; i<ctx->images->n_images; i++ ) {
-		image_feature_list_free(ctx->images->images[i].rflist);
-		ctx->images->images[i].rflist = NULL;
+
+		ImageRecord *image;
+		int j;
+
+		image = &ctx->images->images[i];
+		if ( image->rflist != NULL ) {
+			image_feature_list_free(image->rflist);
+			image->rflist = NULL;
+		}
+
+		for ( j=0; j<image->features->n_features; j++ ) {
+			image->features->features[i].partner = NULL;
+		}
+
 	}
 
 }
@@ -235,7 +235,65 @@ void reproject_cell_to_lattice(ControlContext *ctx) {
 /* Notify that ctx->cell has changed (also need to call displaywindow_update) */
 void reproject_lattice_changed(ControlContext *ctx) {
 
+	int i;
+	int total_reprojected = 0;
+	int total_measured = 0;
+	int partnered_reprojected = 0;
+	int partnered_measured = 0;
+
 	reproject_cell_to_lattice(ctx);
-	displaywindow_update_imagestack(ctx->dw);
+	printf("%i reflections\n", ctx->cell_lattice->n_reflections);
 
+	for ( i=0; i<ctx->images->n_images; i++ ) {
+
+		ImageRecord *image;
+		int j;
+
+		image = &ctx->images->images[i];
+
+		/* Perform relrod projection */
+		image->rflist = reproject_get_reflections(image,
+							ctx->cell_lattice);
+
+		/* Loop over reprojected reflections */
+		for ( j=0; j<image->rflist->n_features; j++ ) {
+
+			double d;
+			ImageFeature f;
+			ImageFeature *p;
+
+			f = image->rflist->features[j];
+			p = image->rflist->features[j].partner;
+
+			if ( p != NULL ) {
+				d = distance(f.x, f.y, p->x, p->y);
+				partnered_reprojected++;
+			}
+		}
+		total_reprojected += image->rflist->n_features;
+
+		/* Loop over measured reflections */
+		for ( j=0; j<image->features->n_features; j++ ) {
+
+			double d;
+			ImageFeature f;
+			ImageFeature *p;
+
+			f = image->features->features[j];
+			p = image->features->features[j].partner;
+
+			if ( p != NULL ) {
+				d = distance(f.x, f.y, p->x, p->y);
+				partnered_measured++;
+			}
+		}
+		total_measured += image->features->n_features;
+	}
+	printf("%i images\n", ctx->images->n_images);
+	printf("%i measured reflections\n", total_measured);
+	printf("%i reprojected reflections\n", total_reprojected);
+	printf("%i partnered measured reflections\n", partnered_measured);
+	printf("%i partnered reprojected reflections\n", partnered_reprojected);
+
+	displaywindow_update_imagestack(ctx->dw);
 }