Enable debugging visualisation

10 files changed, 53 insertions, 41 deletions

diff --git a/src/displaywindow.c b/src/displaywindow.c
index 7abc640..61540b8 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -497,7 +497,7 @@ void displaywindow_update_imagestack(DisplayWindow *dw) {
 		
 		/* Perform relrod projection if necessary */
 		if ( !image->rflist ) {
-			image->rflist = reproject_get_reflections(image, dw->ctx->cell_lattice);
+			image->rflist = reproject_get_reflections(image, dw->ctx->cell_lattice, dw->ctx);
 		}
 		
 		/* Draw the reprojected peaks */
diff --git a/src/intensities.c b/src/intensities.c
index a2e19c3..f1bdf19 100644
--- a/src/intensities.c
+++ b/src/intensities.c
@@ -46,7 +46,7 @@ void intensities_extract(ControlContext *ctx) {
 		ImageRecord *image;
 		
 		image = &ctx->images->images[i];
-		if ( image->rflist == NULL ) image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
+		if ( image->rflist == NULL ) image->rflist = reproject_get_reflections(image, ctx->cell_lattice, ctx);
 		
 		for ( j=0; j<image->rflist->n_features; j++ ) {
 		
diff --git a/src/refine.c b/src/refine.c
index e3bf310..19f9ea3 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -300,7 +300,7 @@ double refine_do_cell(ControlContext *ctx) {
 		
 		if ( !ctx->images->images[i].rflist ) {
 			ctx->images->images[i].rflist = reproject_get_reflections(&ctx->images->images[i],
-										  ctx->cell_lattice);
+										  ctx->cell_lattice, ctx);
 		}
 		
 		for ( j=0; j<ctx->images->images[i].rflist->n_features; j++ ) {
diff --git a/src/refinetest2d.c b/src/refinetest2d.c
index a959bf0..2b09ef4 100644
--- a/src/refinetest2d.c
+++ b/src/refinetest2d.c
@@ -107,6 +107,7 @@ int main(int argc, char *argv[]) {
 	ctx->x_centre = 256;
 	ctx->y_centre = 256;
 	ctx->pixel_size = 5e7;
+	ctx->reflectionlist = reflectionlist_new();
 	image_add(ctx->images, NULL, 512, 512, deg2rad(0.0), ctx);
 	ctx->images->images[0].features = image_feature_list_new();
 	
@@ -121,7 +122,7 @@ int main(int argc, char *argv[]) {
 	cell_real->c.x = 0.0e9;  cell_real->c.y = 0.0e9;  cell_real->c.z = 5.0e9;
 	/* The "real" reflections */
 	reflections_real = reflection_list_from_cell(cell_real);
-	ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real);
+	ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real, ctx);
 //	printf("RT: %i test features generated.\n", ctx->images->images[0].features->n_features);
 	
 	/* The "model" cell to be refined */
@@ -131,13 +132,13 @@ int main(int argc, char *argv[]) {
 	ctx->cell->c.x = 0.1e9;  ctx->cell->c.y = 0.1e9;  ctx->cell->c.z = 5.3e9;
 	ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
 	ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
-	ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice);
+	ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice, ctx);
 	reproject_partner_features(ctx->images->images[0].rflist, &ctx->images->images[0]);
 	
 	refine_do_cell(ctx);
 	image_feature_list_free(ctx->images->images[0].rflist);
 	reflection_list_from_new_cell(ctx->cell_lattice, ctx->cell);
-	ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice);
+	ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice, ctx);
 	
 	fail = check_cell(ctx->cell, cell_real);
 	
diff --git a/src/refinetest3d1.c b/src/refinetest3d1.c
index 4a87fd8..6dfc74a 100644
--- a/src/refinetest3d1.c
+++ b/src/refinetest3d1.c
@@ -114,6 +114,7 @@ int main(int argc, char *argv[]) {
 	ctx->x_centre = 256;
 	ctx->y_centre = 256;
 	ctx->pixel_size = 5e7;
+	ctx->reflectionlist = reflectionlist_new();
 	image_add(ctx->images, NULL, 512, 512, deg2rad(00.0), ctx);
 	ctx->images->images[0].features = image_feature_list_new();
 	image_add(ctx->images, NULL, 512, 512, deg2rad(90.0), ctx);
@@ -133,8 +134,8 @@ int main(int argc, char *argv[]) {
 	cell_real->c.x = 0.0e9;  cell_real->c.y = 0.0e9;  cell_real->c.z = 5.0e9;
 	/* The "real" reflections */
 	reflections_real = reflection_list_from_cell(cell_real);
-	ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real);
-	ctx->images->images[1].features = reproject_get_reflections(&ctx->images->images[1], reflections_real);
+	ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real, ctx);
+	ctx->images->images[1].features = reproject_get_reflections(&ctx->images->images[1], reflections_real, ctx);
 	//ctx->images->images[2].features = reproject_get_reflections(&ctx->images->images[2], reflections_real);
 	
 	/* The "model" cell to be refined */
@@ -143,9 +144,9 @@ int main(int argc, char *argv[]) {
 	ctx->cell->b.x = 0.1e9;  ctx->cell->b.y = 5.0e9;  ctx->cell->b.z = 0.1e9;
 	ctx->cell->c.x = 0.1e9;  ctx->cell->c.y = 0.1e9;  ctx->cell->c.z = 5.0e9;
 	ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
-	ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice);
+	ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice, ctx);
 	reproject_partner_features(ctx->images->images[0].rflist, &ctx->images->images[0]);
-	ctx->images->images[1].rflist = reproject_get_reflections(&ctx->images->images[1], ctx->cell_lattice);
+	ctx->images->images[1].rflist = reproject_get_reflections(&ctx->images->images[1], ctx->cell_lattice, ctx);
 	reproject_partner_features(ctx->images->images[1].rflist, &ctx->images->images[1]);
 	//ctx->images->images[2].rflist = reproject_get_reflections(&ctx->images->images[2], ctx->cell_lattice);
 	//reproject_partner_features(ctx->images->images[2].rflist, &ctx->images->images[2]);
diff --git a/src/refinetest3d2.c b/src/refinetest3d2.c
index 7733106..cffe7f2 100644
--- a/src/refinetest3d2.c
+++ b/src/refinetest3d2.c
@@ -126,6 +126,7 @@ int main(int argc, char *argv[]) {
 	ctx->x_centre = 256;
 	ctx->y_centre = 256;
 	ctx->pixel_size = 5e7;
+	ctx->reflectionlist = reflectionlist_new();
 	image_add(ctx->images, NULL, 512, 512, deg2rad(00.0), ctx);
 	ctx->images->images[0].features = image_feature_list_new();
 	image_add(ctx->images, NULL, 512, 512, deg2rad(90.0), ctx);
@@ -145,9 +146,9 @@ int main(int argc, char *argv[]) {
 	cell_real->c.x = 0.0e9;  cell_real->c.y = 0.0e9;  cell_real->c.z = 5.0e9;
 	/* The "real" reflections */
 	reflections_real = reflection_list_from_cell(cell_real);
-	ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real);
-	ctx->images->images[1].features = reproject_get_reflections(&ctx->images->images[1], reflections_real);
-	ctx->images->images[2].features = reproject_get_reflections(&ctx->images->images[2], reflections_real);
+	ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real, ctx);
+	ctx->images->images[1].features = reproject_get_reflections(&ctx->images->images[1], reflections_real, ctx);
+	ctx->images->images[2].features = reproject_get_reflections(&ctx->images->images[2], reflections_real, ctx);
 	
 	/* The "model" cell to be refined */
 	ctx->cell = malloc(sizeof(Basis));
@@ -155,11 +156,11 @@ int main(int argc, char *argv[]) {
 	ctx->cell->b.x = 0.1e9;  ctx->cell->b.y = 5.0e9;  ctx->cell->b.z = 0.1e9;
 	ctx->cell->c.x = 0.1e9;  ctx->cell->c.y = 0.1e9;  ctx->cell->c.z = 5.0e9;
 	ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
-	ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice);
+	ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice, ctx);
 	reproject_partner_features(ctx->images->images[0].rflist, &ctx->images->images[0]);
-	ctx->images->images[1].rflist = reproject_get_reflections(&ctx->images->images[1], ctx->cell_lattice);
+	ctx->images->images[1].rflist = reproject_get_reflections(&ctx->images->images[1], ctx->cell_lattice, ctx);
 	reproject_partner_features(ctx->images->images[1].rflist, &ctx->images->images[1]);
-	ctx->images->images[2].rflist = reproject_get_reflections(&ctx->images->images[2], ctx->cell_lattice);
+	ctx->images->images[2].rflist = reproject_get_reflections(&ctx->images->images[2], ctx->cell_lattice, ctx);
 	reproject_partner_features(ctx->images->images[2].rflist, &ctx->images->images[2]);
 	
 	refine_do_cell(ctx);
diff --git a/src/refinetest3d3.c b/src/refinetest3d3.c
index 5850b09..65705a8 100644
--- a/src/refinetest3d3.c
+++ b/src/refinetest3d3.c
@@ -127,6 +127,7 @@ int main(int argc, char *argv[]) {
 	ctx->x_centre = 256;
 	ctx->y_centre = 256;
 	ctx->pixel_size = 5e7;
+	ctx->reflectionlist = reflectionlist_new();
 	for ( i=0; i<=90; i++ ) {
 		image_add(ctx->images, NULL, 512, 512, deg2rad(i), ctx);
 		ctx->images->images[i].features = image_feature_list_new();
@@ -144,7 +145,8 @@ int main(int argc, char *argv[]) {
 	/* The "real" reflections */
 	reflections_real = reflection_list_from_cell(cell_real);
 	for ( i=0; i<=90; i++ ) {
-		ctx->images->images[i].features = reproject_get_reflections(&ctx->images->images[i], reflections_real);
+		ctx->images->images[i].features = reproject_get_reflections(&ctx->images->images[i], reflections_real,
+									    ctx);
 	}
 	
 	/* The "model" cell to be refined */
@@ -154,7 +156,8 @@ int main(int argc, char *argv[]) {
 	ctx->cell->c.x = 0.1e9;  ctx->cell->c.y = 0.1e9;  ctx->cell->c.z = 5.0e9;
 	ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
 	for ( i=0; i<ctx->images->n_images; i++ ) {
-		ctx->images->images[i].rflist = reproject_get_reflections(&ctx->images->images[i], ctx->cell_lattice);
+		ctx->images->images[i].rflist = reproject_get_reflections(&ctx->images->images[i], ctx->cell_lattice,
+									  ctx);
 		reproject_partner_features(ctx->images->images[i].rflist, &ctx->images->images[i]);
 	}
 	
diff --git a/src/reflections.h b/src/reflections.h
index 8fafe4b..e51c530 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -22,7 +22,8 @@ typedef enum {
 	REFLECTION_MARKER,
 	REFLECTION_VECTOR_MARKER_1,
 	REFLECTION_VECTOR_MARKER_2,
-	REFLECTION_VECTOR_MARKER_3
+	REFLECTION_VECTOR_MARKER_3,
+	REFLECTION_VECTOR_MARKER_4
 } ReflectionType;
 
 typedef struct reflection_struct {
diff --git a/src/reproject.c b/src/reproject.c
index e1ce13e..3bfc20a 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -45,7 +45,7 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) {
 
 }
 
-ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist) {
+ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist, ControlContext *ctx) {
 
 	ImageFeatureList *flist;
 	Reflection *reflection;
@@ -54,7 +54,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 	double nx, ny, nz, nxt, nyt, nzt;	/* "normal" vector (and calculation intermediates) */
 	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;
+	int done_debug = 0;
+	int i = 0;
 			
 	flist = image_feature_list_new();
 	
@@ -70,7 +71,7 @@ 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);
+	reflection_add(ctx->reflectionlist, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1);
 	
 	/* Determine where "up" is
 	 *	See above. */
@@ -78,7 +79,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 	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);
+	reflection_add(ctx->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2);
 	
 	reflection = reflectionlist->reflections;
 	do {
@@ -125,8 +126,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 			
 				double psi, disc;
 				
-				//reflection_add(ctx->reflectionlist, xl, yl, zl, 1, REFLECTION_GENERATED);
-				//reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_MARKER);
+				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 */
@@ -135,33 +136,36 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 				/* 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 ( (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");
+					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");
+					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");
+					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");
+					if ( (i==20) && !done_debug ) printf("case 4\n");
 				}
 				
-			//	if ( (i==20) && !done_debug ) printf("final psi=%f clockwise from 'up'\n", rad2deg(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));
+				if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +x\n",
+									rad2deg(psi));
 				
 				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;
+				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 */
 				if ( image->fmode == FORMULATION_CLEN ) {
@@ -197,7 +201,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 			} /* else this is the central beam so don't worry about it */
 		
 		}
-		
+		i++;
 		reflection = reflection->next;
 	
 	} while ( reflection );
diff --git a/src/reproject.h b/src/reproject.h
index 93f010a..37eb84d 100644
--- a/src/reproject.h
+++ b/src/reproject.h
@@ -19,7 +19,8 @@
 #include "control.h"
 #include "image.h"
 
-extern ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist);
+extern ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist,
+						   ControlContext *ctx);
 extern void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image);
 extern void reproject_cell_to_lattice(ControlContext *ctx);
 extern void reproject_lattice_changed(ControlContext *ctx);