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

1 files changed, 27 insertions, 14 deletions

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);
 		
 	}