Attempt to make itrans-stat quantitative

Size of ImageDisplayMark circles depend on intensity of reflection


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

5 files changed, 57 insertions, 25 deletions

diff --git a/src/imagedisplay.c b/src/imagedisplay.c
index 1caaf9e..9c3f6cd 100644
--- a/src/imagedisplay.c
+++ b/src/imagedisplay.c
@@ -173,6 +173,7 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 	
 	double scale, xoffs, yoffs;
 	ImageDisplayMark *cur;
+	double max;
 	
 	xoffs = ((double)imagedisplay->drawingarea_width - imagedisplay->view_width) / 2;
 	yoffs = ((double)imagedisplay->drawingarea_height - imagedisplay->view_height) / 2;
@@ -211,6 +212,15 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 			imagedisplay->imagerecord.x_centre * scale,
 			imagedisplay->imagerecord.y_centre * scale + 10);
 	}
+
+	/* Find the maximum intensity */
+	cur = imagedisplay->marks;
+	max = 0.0;
+	while ( cur ) {
+		if ( cur->weight > max ) max = cur->weight;
+		if ( cur->weight < 0.0 ) printf("ID: Warning: ImageDisplayMark with negative weight\n");
+		cur = cur->next;
+	}
 	
 	cur = imagedisplay->marks;
 	while ( cur ) {
@@ -229,10 +239,14 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 		if ( (cur->type == IMAGEDISPLAY_MARK_CIRCLE_1) || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_2)
 			 || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_3) ) {
 		
+			double r;
+		
+			r = 20 * (cur->weight/max);
+			
 			gdk_draw_arc(drawingarea->window, gc, FALSE,
-						xoffs + cur->x*scale - 5,
-						yoffs + imagedisplay->view_height-1-cur->y*scale - 5,
-						11, 11, 0, 64*360);
+						xoffs + cur->x*scale - r,
+						yoffs + imagedisplay->view_height-1-cur->y*scale - r,
+						2*r, 2*r, 0, 64*360);
 			
 		} else if ( (cur->type == IMAGEDISPLAY_MARK_LINE_1) || (cur->type == IMAGEDISPLAY_MARK_LINE_2) ) {
 		
@@ -344,13 +358,14 @@ ImageDisplay *imagedisplay_open(ImageRecord image, const char *title, ImageDispl
 	return imagedisplay_open_with_message(image, title, NULL, flags, NULL, NULL);
 }
 
-void imagedisplay_add_mark(ImageDisplay *imagedisplay, double x, double y, ImageDisplayMarkType type) {
+void imagedisplay_add_mark(ImageDisplay *imagedisplay, double x, double y, ImageDisplayMarkType type, double weight) {
 
 	ImageDisplayMark *new;
 	
 	new = malloc(sizeof(ImageDisplayMark));
 	new->x = x;  new->y = y;
 	new->type = type;
+	new->weight = weight;
 	new->next = NULL;
 	
 	if ( !imagedisplay->marks ) {
@@ -373,6 +388,7 @@ void imagedisplay_add_line(ImageDisplay *imagedisplay, double x1, double y1, dou
 	new->x = x1;  new->y = y1;
 	new->x2 = x2;  new->y2 = y2;
 	new->type = type;
+	new->weight = 1.0;	/* This field makes little sense for a line */
 	new->next = NULL;
 	
 	if ( !imagedisplay->marks ) {
diff --git a/src/imagedisplay.h b/src/imagedisplay.h
index 7b16dc2..681de79 100644
--- a/src/imagedisplay.h
+++ b/src/imagedisplay.h
@@ -44,6 +44,7 @@ typedef struct struct_imagedisplaymark {
 	double				x2;
 	double				y2;
 	ImageDisplayMarkType		type;
+	double				weight;	/* The 'intensity' if this is a peak */
 	
 	struct struct_imagedisplaymark	*next;
 	
@@ -81,7 +82,7 @@ typedef struct imagedisplay_struct {
 extern ImageDisplay *imagedisplay_open(ImageRecord image, const char *title, ImageDisplayFlags flags);
 extern ImageDisplay *imagedisplay_open_with_message(ImageRecord image, const char *title, const char *message,
 						ImageDisplayFlags flags, GCallback mouse_click_func, gpointer callback_data);
-extern void imagedisplay_add_mark(ImageDisplay *imagedisplay, double x, double y, ImageDisplayMarkType type);
+extern void imagedisplay_add_mark(ImageDisplay *imagedisplay, double x, double y, ImageDisplayMarkType type, double weight);
 extern void imagedisplay_add_line(ImageDisplay *imagedisplay, double x1, double y1, double x2, double y2, ImageDisplayMarkType type);
 extern void imagedisplay_force_redraw(ImageDisplay *imagedisplay);
 extern void imagedisplay_put_data(ImageDisplay *imagedisplay, ImageRecord imagerecord);
diff --git a/src/itrans-stat.c b/src/itrans-stat.c
index 67ed96b..fe08a3b 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -329,13 +329,14 @@ static gsl_matrix *itrans_peaksearch_stat_matrix_expand(gsl_matrix *m, int oldsi
 	
 	//printf("me: %d->%d\n",oldsize,newsize);
 	
-	new = gsl_matrix_calloc(2, newsize);
+	new = gsl_matrix_calloc(3, newsize);
 	//printf("me: calloc(2,%d)\n",newsize);
 	for ( j=0; j<oldsize; j++) {
 		if ( j < newsize ) {
 			//printf("me: copying col %d\n",j);
 			gsl_matrix_set(new, 0, j, gsl_matrix_get(m, 0, j));
 			gsl_matrix_set(new, 1, j, gsl_matrix_get(m, 1, j));
+			gsl_matrix_set(new, 2, j, gsl_matrix_get(m, 2, j));
 		}
 	}
 
@@ -353,7 +354,7 @@ static gsl_matrix *itrans_peaksearch_stat_matrix_expand(gsl_matrix *m, int oldsi
  * have to return a pointer to com each time because if the matrix size has to be changed then we need to keep
  * track of the location of the resized instance
  */
-static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double threshold, gsl_matrix *m, gsl_matrix *com, int *com_n, int *com_size) {
+static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double threshold, gsl_matrix *m, gsl_matrix *com, int *com_n, int *com_size, double *val) {
 
 	if ( (i >= 0) && (i < m->size1) ) {
 		if ( (j >= 0) && (j < m->size2) ) {
@@ -361,6 +362,7 @@ static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double
 				if ( gsl_matrix_get(m, i, j) > threshold ) {
 				
 					//printf("dff: found valid point (%d,%d)\n",i,j);
+					*val += gsl_matrix_get(m, i, j);
 					gsl_matrix_set(com, 0, *com_n, i);
 					gsl_matrix_set(com, 1, *com_n, j);
 					*com_n = *com_n + 1;
@@ -369,10 +371,10 @@ static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double
 						*com_size *= 2;
 					}
 					mask[i+j*m->size1] = 1;
-					com = itrans_peaksearch_stat_do_ff(i+1, j, mask,threshold, m, com, com_n, com_size);
-					com = itrans_peaksearch_stat_do_ff(i-1, j, mask,threshold, m, com, com_n, com_size);
-					com = itrans_peaksearch_stat_do_ff(i, j+1, mask,threshold, m, com, com_n, com_size);
-					com = itrans_peaksearch_stat_do_ff(i, j-1, mask,threshold, m, com, com_n, com_size);
+					com = itrans_peaksearch_stat_do_ff(i+1, j, mask,threshold, m, com, com_n, com_size, val);
+					com = itrans_peaksearch_stat_do_ff(i-1, j, mask,threshold, m, com, com_n, com_size, val);
+					com = itrans_peaksearch_stat_do_ff(i, j+1, mask,threshold, m, com, com_n, com_size, val);
+					com = itrans_peaksearch_stat_do_ff(i, j-1, mask,threshold, m, com, com_n, com_size, val);
 					
 				}
 			}
@@ -390,7 +392,8 @@ static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double
  * dependent on implementation (default 4MiB stack for unix?)
  * Implements a crude variable-sized array method which hopefully works
  * Returns a gsl_matrix with x coordinates in row 0 and y
- * coordinates in row 1, which should be of the right length
+ * coordinates in row 1, which should be of the right length.
+ * Intensities (sum of included pixel values) in row 2.
  * Variable count is set to the number of points found
  */
 static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double threshold, int *count) {
@@ -405,7 +408,7 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double thresh
 	mask = calloc(m->size1*m->size2, sizeof(int));
 	size = 32;
 	n = 0;
-	p = gsl_matrix_calloc(2, size);
+	p = gsl_matrix_calloc(3, size);
 
 	//printf("ff: starting loop\n");
 	for ( i=0; i<m->size1; i++ ) {
@@ -413,13 +416,16 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double thresh
 			if ( gsl_matrix_get(m, i, j) > threshold ) {
 				if ( mask[i+j*m->size1] == 0 ) {
 				
+					double val;
+				
 					//printf("ff: found starting point (%d,%d)\n",i,j);
 					com_size = 32;
 					com_n = 0;
 					com_x = com_y = 0.;					
-					com = gsl_matrix_calloc(2, com_size); //this is going to hold the points found for this location
+					com = gsl_matrix_calloc(3, com_size); //this is going to hold the points found for this location
 					//printf("ff: starting floodfill stack\n");
-					com = itrans_peaksearch_stat_do_ff(i, j, mask, threshold, m, com, &com_n, &com_size);
+					val = 0;
+					com = itrans_peaksearch_stat_do_ff(i, j, mask, threshold, m, com, &com_n, &com_size, &val);
 					//printf("ff: ended floodfill stack\n");
 					for ( k=0; k<com_n; k++ ) {
 						com_x += gsl_matrix_get(com, 0, k);
@@ -428,8 +434,11 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double thresh
 					com_x /= com_n;
 					com_y /= com_n;
 					//printf("ff: point CoM (%lf,%lf)\n",com_x,com_y);
+					
+					/* Now add it to the list */
 					gsl_matrix_set(p, 0, n, com_x);
 					gsl_matrix_set(p, 1, n, com_y);
+					gsl_matrix_set(p, 2, n, val);
 					n++;
 					if ( n == size ) {
 						p = itrans_peaksearch_stat_matrix_expand(p, size, size*2);
@@ -486,6 +495,8 @@ static gsl_matrix *itrans_peaksearch_stat_iterate(gsl_matrix *m, unsigned int *c
 		itrans_peaksearch_stat_sigma_filter(m,k);
 		//printf("Iterate: floodfilling\n");
 		p = itrans_peaksearch_stat_floodfill(m, 0, &cur);
+		
+		/* Check for convergence of the number of peaks found */
 		//printf("Iterate: %d points found\n", cur);
 		if ( old < 1.05*cur ) break;
 		old = cur;
@@ -506,7 +517,6 @@ ImageFeatureList *itrans_peaksearch_stat(ImageRecord *imagerecord) {
 	gsl_matrix *m;
 	gsl_matrix *p;
 	int i;
-	double px,py;
 	uint16_t *image = imagerecord->image;
 	ImageFeatureList *flist;
 	
@@ -517,9 +527,12 @@ ImageFeatureList *itrans_peaksearch_stat(ImageRecord *imagerecord) {
 	
 	for ( i=0; i<n_reflections; i++ ) {
 	
-		px = gsl_matrix_get(p,0,i);
-		py = gsl_matrix_get(p,1,i);
-		image_add_feature(flist, px, py, imagerecord, 1.0);
+		double px, py, pi;
+	
+		px = gsl_matrix_get(p, 0, i);
+		py = gsl_matrix_get(p, 1, i);
+		pi = gsl_matrix_get(p, 2, i);
+		image_add_feature(flist, px, py, imagerecord, pi);
 		
 	}
 	
diff --git a/src/refine.c b/src/refine.c
index 6b6eefc..8b6c716 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -172,17 +172,17 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio
 	
 	/* Update the cell */
 	mapping_rotate(dax, day, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("a changed by %f %f %f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9,
+	printf("a changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9,
 										  100*dlx/cell->a.x, 100*dly/cell->a.y, 100*dlz/cell->a.z);
 	cell->a.x += dlx;  cell->a.y += dly;  cell->a.z += dlz;
 	
 	mapping_rotate(dbx, dby, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("b changed by %f %f %f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9,
+	printf("b changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9,
 										  100*dlx/cell->b.x, 100*dly/cell->b.y, 100*dlz/cell->b.z);
 	cell->b.x += dlx;  cell->b.y += dly;  cell->b.z += dlz;
 	
 	mapping_rotate(dcx, dcy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("c changed by %f %f %f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9,
+	printf("c changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9,
 										  100*dlx/cell->c.x, 100*dly/cell->c.y, 100*dlz/cell->c.z);
 	cell->c.x += dlx;  cell->c.y += dly;  cell->c.z += dlz;
 	
@@ -334,7 +334,7 @@ static gint refine_dialog_fit_image(GtkWidget *step_button, ControlContext *ctx)
 		displaywindow_update(ctx->dw);
 		
 		if ( fitted ) {
-			imagedisplay_add_mark(ctx->reproject_id, fitted->x,fitted->y, IMAGEDISPLAY_MARK_CIRCLE_3);
+			imagedisplay_add_mark(ctx->reproject_id, fitted->x,fitted->y, IMAGEDISPLAY_MARK_CIRCLE_3, fitted->intensity);
 		}
 		
 	} else {
diff --git a/src/reproject.c b/src/reproject.c
index 1fd061a..2bb33dc 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -217,13 +217,15 @@ static void reproject_mark_peaks(ControlContext *ctx) {
 		image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
 	}
 	for ( j=0; j<image->rflist->n_features; j++ ) {
-		imagedisplay_add_mark(ctx->reproject_id, image->rflist->features[j].x, image->rflist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_1);
+		imagedisplay_add_mark(ctx->reproject_id, image->rflist->features[j].x, image->rflist->features[j].y,
+					IMAGEDISPLAY_MARK_CIRCLE_1, image->rflist->features[j].intensity);
 	}
 	
 	/* Now draw the original measured peaks */
 	flist = image->features;
 	for ( j=0; j<flist->n_features; j++ ) {
-		imagedisplay_add_mark(ctx->reproject_id, flist->features[j].x, flist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_2);
+		imagedisplay_add_mark(ctx->reproject_id, flist->features[j].x, flist->features[j].y,
+					IMAGEDISPLAY_MARK_CIRCLE_2, image->rflist->features[j].intensity);
 	}
 	
 	/* Now connect partners */