From 5f59200e66b5e8d8b5b6dd963568e1fe6d643b06 Mon Sep 17 00:00:00 2001
From: taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1>
Date: Mon, 27 Aug 2007 22:54:01 +0000
Subject: Working reprojection (NB check directions)

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@76 bf6ca9ba-c028-0410-8290-897cf20841d1
---
 src/reproject.c | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

(limited to 'src/reproject.c')

diff --git a/src/reproject.c b/src/reproject.c
index c6cf2ce..f3c3480 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -15,14 +15,14 @@
 #include "control.h"
 #include "reflections.h"
 
-#define MAX_IMAGE_REFLECTIONS 1024
+#define MAX_IMAGE_REFLECTIONS 8*1024
 
-ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionContext *rctx) {
+ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionContext *rctx, ControlContext *ctx) {
 
 	ImageReflection *refl;
 	Reflection *reflection;
 	int i;
-	double smax = 1e9;
+	double smax = 0.5e9;
 	double tilt, omega;
 	
 	tilt = 2*M_PI*(image.tilt/360);
@@ -60,10 +60,10 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 		c = xl*xl + yl*yl + zl*zl - 2.0*zl/image.lambda;
 		s1 = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
 		s2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a);
-		if ( s1 < s2 ) s = s1; else s = s2;
+		if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2;
 		
 		/* Skip this reflection if s is large */
-		if ( s <= smax ) {
+		if ( fabs(s) <= smax ) {
 		
 			double xddd, yddd, zddd;
 			double xdd, ydd, zdd;
@@ -73,7 +73,9 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 			
 			/* Determine the intersection point */
 			xddd = xl + s*nx;  yddd = yl + s*ny;  zddd = zl + s*nz;
-			printf("intersection at %f,%f,%f\n", xddd, yddd,zddd);
+			reflection_add(ctx->reflectionctx, xl, yl, zl, 1, REFLECTION_CENTRAL);
+			reflection_add(ctx->reflectionctx, xddd, yddd, zddd, 1, REFLECTION_MARKER);
+			
 			/* Invert the image->3D mapping to get the image coordinates */
 			xdd = xddd;
 			ydd = (yddd/cos(tilt) - zddd*tan(tilt)/cos(tilt))/(1+tan(tilt)*tan(tilt));
@@ -85,10 +87,11 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 			
 			if ( image.fmode == FORMULATION_CLEN ) {
 				psi = atan2(-yd, xd);
-				theta = acos(1-zd*image.lambda);
+				theta = acos(1+zd*image.lambda);
 				x = image.camera_len*sin(theta)*cos(psi);
 				y = image.camera_len*sin(theta)*sin(psi);
 				x *= image.resolution;
+				y *= image.resolution;
 			} else if ( image.fmode == FORMULATION_PIXELSIZE ) {
 				x = xd / image.pixel_size;
 				y = yd / image.pixel_size;
@@ -117,7 +120,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 				printf("Reflection %i at %i,%i\n", i, refl[i-1].x, refl[i-1].y);
 			
 			} else {
-				fprintf(stderr, "Reflection failed sanity test (x=%f, y=%f)\n", x, y);
+				//fprintf(stderr, "Reflection failed sanity test (x=%f, y=%f)\n", x, y);
 			}
 		
 		}
-- 
cgit v1.2.3