Debodge reprojection and ImageDisplay

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@123 bf6ca9ba-c028-0410-8290-897cf20841d1

4 files changed, 114 insertions, 59 deletions

diff --git a/src/displaywindow.c b/src/displaywindow.c
index 1fef6ce..a1e2a1a 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -282,6 +282,46 @@ static gint displaywindow_gl_motion_notify(GtkWidget *widget, GdkEventMotion *ev
 			ADD_VERTEX					\
 		}							\
 	}
+	
+#define DRAW_POINTER						\
+	glBegin(GL_LINES);					\
+		glVertex3f(50, 0.0, 0.0);			\
+		glVertex3f(-50, 0.0, 0.0);			\
+	glEnd();						\
+	glBegin(GL_TRIANGLE_FAN);				\
+		glVertex3f(50.0, 0.0, 0.0);			\
+		glNormal3f(1.0, 0.0, 0.0);			\
+		glVertex3f(50.0-5.0, 1.0, 1.0);			\
+		glNormal3f(-2.5, 1.0, 1.0);			\
+		glVertex3f(50.0-5.0, -1.0, 1.0);		\
+		glNormal3f(-2.5, -1.0, 1.0);			\
+		glVertex3f(50.0-5.0, -1.0, -1.0);		\
+		glNormal3f(-2.5, -1.0, -1.0);			\
+		glVertex3f(50.0-5.0, 1.0, -1.0);		\
+		glNormal3f(-2.5, 1.0, -1.0);			\
+		glVertex3f(50.0-5.0, 1.0, 1.0);			\
+		glNormal3f(-2.5, 1.0, 1.0);			\
+	glEnd();
+
+#define DRAW_SHORT_POINTER					\
+	glBegin(GL_LINES);					\
+		glVertex3f(0.0, 0.0, 0.0);			\
+		glVertex3f(10.0, 0.0, 0.0);			\
+	glEnd();						\
+	glBegin(GL_TRIANGLE_FAN);				\
+		glVertex3f(10.0, 0.0, 0.0);			\
+		glNormal3f(1.0, 0.0, 0.0);			\
+		glVertex3f(10.0-2.0, 0.2, 0.2);			\
+		glNormal3f(-2.5, 0.2, 0.2);			\
+		glVertex3f(10.0-2.0, -0.2, 0.2);		\
+		glNormal3f(-2.5, -0.2, 0.2);			\
+		glVertex3f(10.0-2.0, -0.2, -0.2);		\
+		glNormal3f(-2.5, -0.2, -0.2);			\
+		glVertex3f(10.0-2.0, 0.2, -0.2);		\
+		glNormal3f(-2.5, 0.2, -0.2);			\
+		glVertex3f(10.0-2.0, 0.2, 0.2);			\
+		glNormal3f(-2.5, 0.2, 0.2);			\
+	glEnd();
 
 static void displaywindow_gl_create_list(DisplayWindow *dw) {
 
@@ -435,19 +475,20 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 	glEnd();
 	
 	/* x, y, z pointers */
-	glBegin(GL_LINES);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, red);
-		glVertex3f(0.0, 0.0, 0.0);
-		glVertex3f(5.0, 0.0, 0.0);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, green);
-		glVertex3f(0.0, 0.0, 0.0);
-		glVertex3f(0.0, 5.0, 0.0);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, bblue);
-		glVertex3f(0.0, 0.0, 0.0);
-		glVertex3f(0.0, 0.0, 5.0);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
-	glEnd();
+	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
+	DRAW_SHORT_POINTER
+	
+	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, green);
+	glPushMatrix();
+	glRotatef(90.0, 0.0, 0.0, 1.0);
+	DRAW_SHORT_POINTER
+	glPopMatrix();
+	
+	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, bblue);
+	glPushMatrix();
+	glRotatef(-90.0, 0.0, 1.0, 0.0);
+	DRAW_SHORT_POINTER
+	glPopMatrix();
 	
 	/* Tilt axis */
 	glPushMatrix();
@@ -458,30 +499,10 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 	glRotatef(ctx->images[0].omega, 0.0, 0.0, 1.0);
 	glPolygonOffset(1.0, 1.0);
 	glEnable(GL_POLYGON_OFFSET_LINE);
-	glBegin(GL_LINES);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, yellow);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
-		glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
-		glVertex3f(50, 0.0, 0.0);
-		glVertex3f(-50, 0.0, 0.0);
-	glEnd();
-	glBegin(GL_TRIANGLE_FAN);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, yellow);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
-		glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
-		glVertex3f(50.0, 0.0, 0.0);
-		glNormal3f(1.0, 0.0, 0.0);
-		glVertex3f(50.0-5.0, 1.0, 1.0);
-		glNormal3f(-2.5, 1.0, 1.0);
-		glVertex3f(50.0-5.0, -1.0, 1.0);
-		glNormal3f(-2.5, -1.0, 1.0);
-		glVertex3f(50.0-5.0, -1.0, -1.0);
-		glNormal3f(-2.5, -1.0, -1.0);
-		glVertex3f(50.0-5.0, 1.0, -1.0);
-		glNormal3f(-2.5, 1.0, -1.0);
-		glVertex3f(50.0-5.0, 1.0, 1.0);
-		glNormal3f(-2.5, 1.0, 1.0);
-	glEnd();
+	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, yellow);
+	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
+	glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
+	DRAW_POINTER
 	glDisable(GL_POLYGON_OFFSET_LINE);
 	glPopMatrix();
 	
@@ -509,7 +530,7 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 			glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
 			glBegin(GL_LINES);
 			glVertex3f(0.0, 0.0, 0.0);
-			glVertex3f(reflection->x, reflection->y, reflection->z);
+			glVertex3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9);
 			glEnd();
 			
 		}
@@ -521,7 +542,7 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 			glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
 			glBegin(GL_LINES);
 			glVertex3f(0.0, 0.0, 0.0);
-			glVertex3f(reflection->x, reflection->y, reflection->z);
+			glVertex3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9);
 			glEnd();
 			
 		}
@@ -533,7 +554,19 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 			glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
 			glBegin(GL_LINES);
 			glVertex3f(0.0, 0.0, 0.0);
-			glVertex3f(reflection->x, reflection->y, reflection->z);
+			glVertex3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9);
+			glEnd();
+			
+		}
+		
+		if ( reflection->type == REFLECTION_VECTOR_MARKER_4 ) {
+		
+			glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, purple);
+			glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
+			glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
+			glBegin(GL_LINES);
+			glVertex3f(0.0, 0.0, 0.0);
+			glVertex3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9);
 			glEnd();
 			
 		}
diff --git a/src/imagedisplay.c b/src/imagedisplay.c
index cfb0c74..8647f74 100644
--- a/src/imagedisplay.c
+++ b/src/imagedisplay.c
@@ -202,7 +202,7 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 	cur = imagedisplay->marks;
 	while ( cur ) {
 		gdk_draw_arc(drawingarea->window, imagedisplay->gc_marks, FALSE,
-						xoffs + imagedisplay->view_height-cur->x*scale - 5,
+						xoffs + cur->x*scale - 5,
 						yoffs + imagedisplay->view_height-cur->y*scale - 5,
 						11, 11, 0, 64*360);
 		cur = cur->next;
diff --git a/src/reflections.h b/src/reflections.h
index 253a110..2dd363c 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -23,7 +23,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 78cced6..06611cd 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -23,11 +23,12 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	ImageReflection *refl;
 	Reflection *reflection;
 	int i;
-	double smax = 0.5e9;
+	double smax = 0.25e9;
 	double tilt, omega;
 	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;
 			
 	tilt = deg2rad(image.tilt);
 	omega = deg2rad(image.omega);
@@ -38,7 +39,8 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	
 	reflection = rctx->reflections;
 	
-	/* Calculate the (normalised) incident electron wavevector */
+	/* Calculate the (normalised) incident electron wavevector
+	 *	n is rotated "into" the reconstruction, so only one omega step. */
 	nxt = 0.0;  nyt = 0.0;  nzt = 1.0;
 	nx = nxt;  ny = cos(tilt)*nyt + sin(tilt)*nzt;  nz = -sin(tilt)*nyt + cos(tilt)*nzt;
 	nxt = nx;  nyt = ny;  nzt = nz;
@@ -48,14 +50,13 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	kz = nz / image.lambda;	/* This is the centre of the Ewald sphere */
 	reflection_add(ctx->reflectionctx, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1);
 	
-	/* Determine where "up" is */
+	/* 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->reflectionctx, ux*50, uy*50, uz*50, 1, REFLECTION_VECTOR_MARKER_2);
-	
-	printf("angle between k and u is %f\n", rad2deg(angle_between(kx, ky, kz, ux, uy, uz)));
+	reflection_add(ctx->reflectionctx, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2);
 	
 	do {
 	
@@ -106,29 +107,49 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 				reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_MARKER);
 				
 				/* Calculate azimuth of point in image (clockwise from "up", will be changed later) */
-				cx = nz*yi - ny*zi;  cy = nx*zi - nz*xi;  cz = ny*xi - nx*yi;  /* c = (-n) x i */
+				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 = nz*cy - ny*cz;  ry = nx*cz - nz*cx;  rz = ny*cx - nx*cy;  /* r = (-n) x [(-n) x i] */
+				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 ( (psi >= M_PI_2) && (disc >= M_PI_2) ) psi -= M_PI_2;		/* Case 1 */
-				else if ( (psi >= M_PI_2) && (disc < M_PI_2) ) psi = 3*M_PI_2 - psi;	/* Case 2 */
-				else if ( (psi < M_PI_2) && (disc < M_PI_2) ) psi = 3*M_PI_2 - psi;	/* Case 3 */
-				else if ( (psi < M_PI_2) && (disc >= M_PI_2) ) psi = 3*M_PI_2 + psi;	/* Case 4 */
+			//	if ( (i==20) && !done_debug ) {
+			//		reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_VECTOR_MARKER_3);
+			//		reflection_add(ctx->reflectionctx, cx, cy, cz, 1, REFLECTION_VECTOR_MARKER_4);
+			//		reflection_add(ctx->reflectionctx, 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");
+				}
+				
+			//	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));
 				
 				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 */
 				if ( image.fmode == FORMULATION_CLEN ) {
-					x = image.camera_len*sin(theta)*cos(psi);
-					y = image.camera_len*sin(theta)*sin(psi);
+					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 ) {
-					x = sin(theta)*cos(psi) / image.lambda;
-					y = sin(theta)*sin(psi) / image.lambda;
+					x = tan(theta)*cos(psi) / image.lambda;
+					y = tan(theta)*sin(psi) / image.lambda;
 					x /= image.pixel_size;
 					y /= image.pixel_size;
 				} else {
@@ -139,7 +160,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 				/* Adjust centre */
 				x += image.x_centre;
 				y += image.y_centre;
-				x = image.width - 1 - x;	/* Debodge me */
+				
 				/* Sanity check */
 				if ( (x>=0) && (x<image.width) && (y>=0) && (y<image.height) ) {