Refactor image handling code

Remove itrans-lsq


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

1 files changed, 45 insertions, 55 deletions

diff --git a/src/reproject.c b/src/reproject.c
index 15a3ca9..3604481 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -17,14 +17,12 @@
 #include "utils.h"
 #include "imagedisplay.h"
 #include "displaywindow.h"
+#include "image.h"
 
-#define MAX_IMAGE_REFLECTIONS 8*1024
+ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist, ControlContext *ctx) {
 
-ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionList *reflectionlist, ControlContext *ctx) {
-
-	ImageReflection *refl;
+	ImageFeatureList *flist;
 	Reflection *reflection;
-	int i;
 	double smax = 0.1e9;
 	double tilt, omega;
 	double nx, ny, nz, nxt, nyt, nzt;	/* "normal" vector (and calculation intermediates) */
@@ -32,8 +30,10 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	double ux, uy, uz, uxt, uyt, uzt;	/* "up" vector (and calculation intermediates) */
 	//int done_debug = 0;
 			
-	tilt = deg2rad(image.tilt);
-	omega = deg2rad(image.omega);
+	flist = image_feature_list_new();
+	
+	tilt = deg2rad(image->tilt);
+	omega = deg2rad(image->omega);
 	
 	/* Calculate the (normalised) incident electron wavevector
 	 *	n is rotated "into" the reconstruction, so only one omega step. */
@@ -41,9 +41,9 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	nx = nxt;  ny = cos(tilt)*nyt + sin(tilt)*nzt;  nz = -sin(tilt)*nyt + cos(tilt)*nzt;
 	nxt = nx;  nyt = ny;  nzt = nz;
 	nx = nxt*cos(-omega) + nyt*sin(-omega);  ny = -nxt*sin(-omega) + nyt*cos(-omega);  nz = nzt;
-	kx = nx / image.lambda;
-	ky = ny / image.lambda;
-	kz = nz / image.lambda;	/* This is the centre of the Ewald sphere */
+	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
@@ -54,9 +54,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	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);
 	
-	refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection));
 	reflection = reflectionlist->reflections;
-	i = 0;
 	do {
 	
 		double xl, yl, zl;
@@ -71,11 +69,11 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 		/* Next, solve the relrod equation to calculate the excitation error */
 		a = 1.0;
 		b = -2.0*(xl*nx + yl*ny + zl*nz);
-		c = xl*xl + yl*yl + zl*zl - 1.0/(image.lambda*image.lambda);
+		c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda);
 		A1 = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
 		A2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a);
-		s1 = 1.0/image.lambda - A1;
-		s2 = 1.0/image.lambda - A2;
+		s1 = 1.0/image->lambda - A1;
+		s2 = 1.0/image->lambda - A2;
 		if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2;
 		
 		/* Skip this reflection if s is large */
@@ -140,39 +138,31 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 			//	if ( (i==20) && !done_debug ) done_debug = 1;
 				
 				/* Calculate 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 ) {
-					x = tan(theta)*cos(psi) / image.lambda;
-					y = tan(theta)*sin(psi) / image.lambda;
-					x /= image.pixel_size;
-					y /= image.pixel_size;
+				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 ) {
+					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");
 					return NULL;
 				}
 				
 				/* Adjust centre */
-				x += image.x_centre;
-				y += image.y_centre;
+				x += image->x_centre;
+				y += image->y_centre;
 				
 				/* Sanity check */
-				if ( (x>=0) && (x<image.width) && (y>=0) && (y<image.height) ) {
+				if ( (x>=0) && (x<image->width) && (y>=0) && (y<image->height) ) {
 				
 					/* Record the reflection */
-					refl[i].x = x;
-					refl[i].y = y;
-					i++;
-					
-					if ( i > MAX_IMAGE_REFLECTIONS ) {
-						fprintf(stderr, "Too many reflections\n");
-						break;
-					}
-					
-					//printf("Reflection %i at %i,%i\n", i, refl[i-1].x, refl[i-1].y);
+					image_add_feature(flist, x, y, image, reflection->intensity);
+					//printf("Reflection at %f, %f\n", x, y);
 				
 				} /* else it's outside the picture somewhere */
 				
@@ -184,30 +174,28 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 
 	} while ( reflection );
 	
-	//printf("Found %i reflections in image\n", i);
-	*n = i;
-	return refl;
+	return flist;
 
 }
 
 static gint reproject_clicked(GtkWidget *widget, GdkEventButton *event, ControlContext *ctx) {
 
-	ImageReflection *refl;
-	size_t n, j;
+	ImageFeatureList *flist;
+	size_t j;
 	
 	ctx->reproject_cur_image++;
-	if ( ctx->reproject_cur_image == ctx->n_images ) ctx->reproject_cur_image = 0;
+	if ( ctx->reproject_cur_image == ctx->images->n_images ) ctx->reproject_cur_image = 0;
 	
 	imagedisplay_clear_circles(ctx->reproject_id);
 	reflectionlist_clear_markers(ctx->reflectionlist);
 	
-	refl = reproject_get_reflections(ctx->images[ctx->reproject_cur_image], &n, ctx->cell_lattice, ctx);
-	for ( j=0; j<n; j++ ) {
-		imagedisplay_mark_circle(ctx->reproject_id, refl[j].x, refl[j].y);
+	flist = reproject_get_reflections(&ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice, ctx);
+	for ( j=0; j<flist->n_features; j++ ) {
+		imagedisplay_mark_circle(ctx->reproject_id, flist->features[j].x, flist->features[j].y);
 	}
+	image_feature_list_free(flist);
 	
-	imagedisplay_put_data(ctx->reproject_id, ctx->images[ctx->reproject_cur_image]);
-	displaywindow_update(ctx->dw);
+	imagedisplay_put_data(ctx->reproject_id, ctx->images->images[ctx->reproject_cur_image]);
 	
 	return 0;
 
@@ -215,8 +203,8 @@ static gint reproject_clicked(GtkWidget *widget, GdkEventButton *event, ControlC
 
 void reproject_open(ControlContext *ctx) {
 	
-	size_t n, j;
-	ImageReflection *refl;
+	size_t j;
+	ImageFeatureList *flist;
 	
 	if ( !ctx->cell ) {
 		printf("RP: No current cell\n");
@@ -233,13 +221,15 @@ void reproject_open(ControlContext *ctx) {
 	ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
 	ctx->reproject_cur_image = 0;
 	
-	ctx->reproject_id = imagedisplay_open_with_message(ctx->images[ctx->reproject_cur_image], "Current Image", "Click to change image",
+	ctx->reproject_id = imagedisplay_open_with_message(ctx->images->images[ctx->reproject_cur_image], "Current Image", "Click to change image",
 				IMAGEDISPLAY_SHOW_CENTRE | IMAGEDISPLAY_SHOW_TILT_AXIS, G_CALLBACK(reproject_clicked), ctx);
 	
-	refl = reproject_get_reflections(ctx->images[ctx->reproject_cur_image], &n, ctx->cell_lattice, ctx);
-	for ( j=0; j<n; j++ ) {
-		imagedisplay_mark_circle(ctx->reproject_id, refl[j].x, refl[j].y);
+	flist = reproject_get_reflections(&ctx->images->images[ctx->reproject_cur_image], ctx->cell_lattice, ctx);
+	for ( j=0; j<flist->n_features; j++ ) {
+		imagedisplay_mark_circle(ctx->reproject_id, flist->features[j].x, flist->features[j].y);
 	}
+	printf("%i reflections\n", flist->n_features);
+	image_feature_list_free(flist);
 	
 }