Line wrapping, whitespace, remove debug

1 files changed, 150 insertions, 159 deletions

diff --git a/src/itrans-stat.c b/src/itrans-stat.c
index fe08a3b..ee41817 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -25,24 +25,25 @@
 /* Renormalise a gsl_matrix to 0->1
  * Re-written to use gsl_matrix native functions */
 static void matrix_renormalise(gsl_matrix *m) {
-	
+
 	double max,min;
-	
+
 	gsl_matrix_minmax(m,&min,&max);
 	gsl_matrix_add_constant(m,0.-min);
 	gsl_matrix_scale(m,1./(max-min));
 
 }
 
-/* 
+/*
  * Create a gsl_matrix for performing image operations on
  * from a raw image and control context
  * Converting to gsl_matrix because many of the required operations
  * can be done as matrices to save effort
  * Renormalises matrix to 0->1
  */
-static gsl_matrix *itrans_peaksearch_stat_createImageMatrix(uint16_t *image, int width, int height) {
-
+static gsl_matrix *itrans_peaksearch_stat_createImageMatrix(uint16_t *image,
+							int width, int height)
+{
 	gsl_matrix *raw;
 	int i, j;
 
@@ -53,14 +54,13 @@ static gsl_matrix *itrans_peaksearch_stat_createImageMatrix(uint16_t *image, int
 		}
 	}
 	matrix_renormalise(raw);
-	
-	return raw;
 
+	return raw;
 }
 
 /* Find and return the mean value of the matrix */
-static double itrans_peaksearch_stat_image_mean(gsl_matrix *m) {
-
+static double itrans_peaksearch_stat_image_mean(gsl_matrix *m)
+{
 	int i, j;
 	double mean = 0.;
 
@@ -71,34 +71,33 @@ static double itrans_peaksearch_stat_image_mean(gsl_matrix *m) {
 	}
 
 	return mean / (m->size1 * m->size2);
-
 }
 
 /*
  * Return the standard deviation, sigma, of the matrix values
  * \sqrt(\sum((x-mean)^2)/n)
  */
-static double itrans_peaksearch_stat_image_sigma(gsl_matrix *m, double mean) {
-
+static double itrans_peaksearch_stat_image_sigma(gsl_matrix *m, double mean)
+{
 	int i,j;
 	double diff2 = 0;
 
 	for ( i=0; i<m->size1; i++ ) {
 		for ( j=0; j<m->size2; j++ ) {
-			diff2 += (gsl_matrix_get(m,i,j)-mean)*(gsl_matrix_get(m,i,j)-mean);
+			diff2 += (gsl_matrix_get(m,i,j)-mean)
+				 * (gsl_matrix_get(m,i,j)-mean);
 		}
 	}
-	
-	return sqrt(diff2/(m->size1 * m->size2));
 
+	return sqrt(diff2/(m->size1 * m->size2));
 }
 
 /*
  * Filter all pixels with value < mean + k*sigma to 0
  * Set all matching pixels to 1 d
  */
-static void itrans_peaksearch_stat_sigma_filter(gsl_matrix *m, double k) {
-
+static void itrans_peaksearch_stat_sigma_filter(gsl_matrix *m, double k)
+{
 	double mean,sigma;
 	int i,j;
 
@@ -114,7 +113,6 @@ static void itrans_peaksearch_stat_sigma_filter(gsl_matrix *m, double k) {
 			}
 		}
 	}
-	
 }
 
 /*
@@ -122,26 +120,25 @@ static void itrans_peaksearch_stat_sigma_filter(gsl_matrix *m, double k) {
  *
  * TODO: Use a mask instead of calculating valid points
  */
-static double itrans_peaksearch_stat_circle_mean(gsl_matrix *m, int x, int y, int r, gsl_matrix *mask) {
-
+static double itrans_peaksearch_stat_circle_mean(gsl_matrix *m, int x, int y,
+						 int r, gsl_matrix *mask)
+{
 	double mean = 0.;
 	int i, j;
 	int n = 0;
 
 	for ( i=x-r; i<=x+r; i++ ) {
 		for ( j=y-r; j<=y+r; j++ ) {
-			//printf("cm: ij=(%d,%d) mask=(%d,%d)\n",i,j,i-x+r,j-y+r);
 			if ( gsl_matrix_get(mask,i-x+r,j-y+r)>0. ) {
 				mean += gsl_matrix_get(m, i, j);
-				//printf("cm: (%d,%d) mean=%lf val=%lf\n",i,j,mean,gsl_matrix_get(m,i,j));
 				n++;
 			}
 		}
 	}
-	
+
 	//printf("cm: (%d,%d) summean=%lf n=%d\n",x,y,mean,n);
 	return mean/n;
-	
+
 }
 
 /*
@@ -149,8 +146,9 @@ static double itrans_peaksearch_stat_circle_mean(gsl_matrix *m, int x, int y, in
  *
  * TODO: Use a mask instead of calculating valid points
  */
-static double itrans_peaksearch_stat_circle_sigma(gsl_matrix *m, int x, int y, int r, gsl_matrix *mask, double mean) {
-
+static double itrans_peaksearch_stat_circle_sigma(gsl_matrix *m, int x, int y,
+					int r, gsl_matrix *mask, double mean)
+{
 	double diff2 = 0.;
 	int i, j;
 	int n = 0;
@@ -158,21 +156,21 @@ static double itrans_peaksearch_stat_circle_sigma(gsl_matrix *m, int x, int y, i
 	for ( i=x-r; i<=x+r; i++ ) {
 		for ( j=y-r; j<=y+r; j++ ) {
 			if ( gsl_matrix_get(mask, i-x+r, j-y+r) > 0 ) {
-				diff2 += (gsl_matrix_get(m,i,j)-mean)*(gsl_matrix_get(m,i,j)-mean);
+				diff2 += (gsl_matrix_get(m,i,j)-mean)
+					 * (gsl_matrix_get(m,i,j)-mean);
 				n++;
 			}
 		}
 	}
-	
-	return sqrt(diff2/n);
 
+	return sqrt(diff2/n);
 }
 
 /*
  * Calculate a circular mask to save recalculating it every time
  */
-static gsl_matrix *itrans_peaksearch_stat_circle_mask(int r) {
-
+static gsl_matrix *itrans_peaksearch_stat_circle_mask(int r)
+{
 	gsl_matrix *m;
 	int i,j;
 
@@ -182,12 +180,11 @@ static gsl_matrix *itrans_peaksearch_stat_circle_mask(int r) {
 			if ( sqrt((r-i)*(r-i)+(r-j)*(r-j)) <= r ) {
 				 gsl_matrix_set(m, i, j, 1.);
 			}
-			
+
 		}
 	}
-	
+
 	return m;
-	
 }
 
 /*
@@ -206,59 +203,60 @@ static gsl_matrix *itrans_peaksearch_stat_circle_mask(int r) {
  *
  * TODO: Pass a mask to the ancilliary functions instead of having
  * them calculate several hundred million sqrts
- * 
- * self-referencing problem being dealt with - output being written onto the array before the next point it computed
- * problem carried over from the OO version where a new object was created by each operation
+ *
+ * self-referencing problem being dealt with - output being written onto the
+ * array before the next point it computed problem carried over from the OO
+ * version where a new object was created by each operation
  */
-static void itrans_peaksearch_stat_local_sigma_filter(gsl_matrix *m, int r, double k) {
-
+static void itrans_peaksearch_stat_local_sigma_filter(gsl_matrix *m, int r,
+							double k)
+{
 	double mean,sigma;
 	double local;
 	gsl_matrix *mask;
 	gsl_matrix *new;
 	int i,j;
-	//int interval;
 
 	mask = itrans_peaksearch_stat_circle_mask(r);
 	new = gsl_matrix_alloc(m->size1, m->size2);
 
-	//interval = (m->size1-r)/20;
-	//printf("lsf: starting loop\n");
-	//printf("lsf: ");
-	//for (i=r;i<m->size1-r;i++) {
-	//	for (j=r;j<m->size2-r;j++) {
-
 	for ( i=0; i<m->size1; i++ ) {
 		for ( j=0; j<m->size2; j++ ) {
-			if ( ((i >= r) && (i < m->size1-r)) && ((j >= r) && (j < m->size2-r)) ) {
-			
-				//printf("lsf: evaluating (%d,%d)\n",i,j);
-				mean = itrans_peaksearch_stat_circle_mean(m, i, j, r, mask);
-				//printf("lsf: mean=%lf",mean);
-				sigma = itrans_peaksearch_stat_circle_sigma(m, i, j, r, mask, mean);
-				//printf(" sigma=%lf",sigma);
+
+			if ( ((i >= r) && (i < m->size1-r))
+			      && ((j >= r) && (j < m->size2-r)) ) {
+
+				mean = itrans_peaksearch_stat_circle_mean(m,
+							i, j, r, mask);
+				sigma = itrans_peaksearch_stat_circle_sigma(m,
+							i, j, r, mask, mean);
+
 				local = gsl_matrix_get(m, i, j);
-				local += gsl_matrix_get(m, i+1, j) + gsl_matrix_get(m, i-1, j) + gsl_matrix_get(m, i, j+1) + gsl_matrix_get(m, i, j-1);
-				local += .5*(gsl_matrix_get(m, i+1, j+1) + gsl_matrix_get(m, i-1, j+1) + gsl_matrix_get(m, i+1, j-1) + gsl_matrix_get(m, i-1, j-1));
+				local += gsl_matrix_get(m, i+1, j)
+					+ gsl_matrix_get(m, i-1, j)
+					+ gsl_matrix_get(m, i, j+1)
+					+ gsl_matrix_get(m, i, j-1);
+				local += .5*(gsl_matrix_get(m, i+1, j+1)
+					+ gsl_matrix_get(m, i-1, j+1)
+					+ gsl_matrix_get(m, i+1, j-1)
+					+ gsl_matrix_get(m, i-1, j-1));
 				local /= 7.;
-				//printf(" local=%lf\n",local);
+
 				if ( local > mean+k*sigma ) {
 					gsl_matrix_set(new, i, j, 1.);
 				} else {
 					gsl_matrix_set(new, i, j, 0.);
 				}
-				
+
 			} else {
 				gsl_matrix_set(new, i, j, 0.);
 			}
+
 		}
-		//if (i % interval == 0) printf(".");
 	}
-	
-	//printf("done\n");
+
 	gsl_matrix_memcpy(m, new);
 	gsl_matrix_free(new);
-
 }
 
 /*
@@ -268,11 +266,12 @@ static void itrans_peaksearch_stat_local_sigma_filter(gsl_matrix *m, int r, doub
  * to be square and should probably be normalised.
  * Don't do daft things like rectangular kernels or kernels larger
  * than the image - this doesn't check such things and will cry.
- * 
+ *
  * Also suffers from self-reference problem
  */
-static void itrans_peaksearch_stat_apply_kernel(gsl_matrix *m, gsl_matrix *kernel) {
-
+static void itrans_peaksearch_stat_apply_kernel(gsl_matrix *m,
+						gsl_matrix *kernel)
+{
 	int size;
 	int half;
 	gsl_matrix *l;
@@ -287,100 +286,102 @@ static void itrans_peaksearch_stat_apply_kernel(gsl_matrix *m, gsl_matrix *kerne
 
 	for ( i=0; i<m->size1; i++ ) {
 		for ( j=0; j<m->size2; j++ ) {
-			if ( ((i >= half) && (i < m->size1-half)) && ((j >= half) && (j < m->size2-half)) ) {
-			
-				lv = gsl_matrix_submatrix(m, i-half, j-half, size, size);
+			if ( ((i >= half) && (i < m->size1-half))
+				&& ((j >= half) && (j < m->size2-half)) ) {
+
+				lv = gsl_matrix_submatrix(m, i-half, j-half,
+								size, size);
 				l = &lv.matrix;
 				val = 0.;
 				for ( x=0; x<size; x++ ) {
 					for ( y=0; y<size; y++ ) {
-						val += gsl_matrix_get(l, x, y)*gsl_matrix_get(kernel, x, y);
+						val += gsl_matrix_get(l, x, y)
+						 * gsl_matrix_get(kernel, x, y);
 					}
 				}
-				//gsl_matrix_free(l);
 				gsl_matrix_set(new,i,j,val);
-				
+
 			}
 		}
 	}
-	
+
 	gsl_matrix_memcpy(m,new);
 	gsl_matrix_free(new);
-	
 }
 
 /* Generate the simplist possible kernel - a flat one */
-static gsl_matrix *itrans_peaksearch_stat_generate_flat_kernel(int half) {
-
+static gsl_matrix *itrans_peaksearch_stat_generate_flat_kernel(int half)
+{
 	gsl_matrix *k;
 
 	k = gsl_matrix_alloc(2*half+1,2*half+1);
 	gsl_matrix_set_all(k,1./((2*half+1)*(2*half+1)));
 
 	return k;
-
 }
 
 /* Expands or contracts a gsl_matrix by copying the columns to a new one */
-static gsl_matrix *itrans_peaksearch_stat_matrix_expand(gsl_matrix *m, int oldsize, int newsize) {
-
+static gsl_matrix *itrans_peaksearch_stat_matrix_expand(gsl_matrix *m,
+				unsigned int oldsize, unsigned int newsize)
+{
 	gsl_matrix *new;
 	int j;
-	
-	//printf("me: %d->%d\n",oldsize,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));
 		}
 	}
 
-	//printf("me: freeing old matrix\n");
 	gsl_matrix_free(m);
 
-	//printf("me: new s1=%d s2=%d\n",new->size1,new->size2);
-	//printf("me: m s1=%d s2=%d\n",m->size1,m->size2);
 	return new;
-	
 }
 
 /*
- * Stack-based flood-fill iteration routine 
- * 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
+ * Stack-based flood-fill iteration routine
+ * 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, double *val) {
-
+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) ) {
 			if ( mask[i+j*m->size1] == 0 ) {
 				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;
 					if ( *com_n == *com_size ) {
-						com = itrans_peaksearch_stat_matrix_expand(com, *com_size, *com_size*2);
+
+	com = itrans_peaksearch_stat_matrix_expand(com, *com_size, *com_size*2);
+
 						*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, 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);
-					
+
+	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);
+
 				}
 			}
 		}
 	}
-	
+
 	return com;
 
 }
@@ -396,8 +397,9 @@ static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double
  * 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) {
-
+static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m,
+					double threshold, int *count)
+{
 	int *mask;
 	int size, com_size, i, j, k, n;
 	int com_n;
@@ -410,41 +412,40 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double thresh
 	n = 0;
 	p = gsl_matrix_calloc(3, size);
 
-	//printf("ff: starting loop\n");
 	for ( i=0; i<m->size1; i++ ) {
 		for ( j=0; j<m->size2; j++ ) {
 			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(3, com_size); //this is going to hold the points found for this location
-					//printf("ff: starting floodfill stack\n");
+					com_x = com_y = 0.;
+					com = gsl_matrix_calloc(3, 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");
+		com = itrans_peaksearch_stat_do_ff(i, j, mask, threshold, m,
+						com, &com_n, &com_size, &val);
 					for ( k=0; k<com_n; k++ ) {
-						com_x += gsl_matrix_get(com, 0, k);
-						com_y += gsl_matrix_get(com, 1, k);
+						com_x += gsl_matrix_get(com,
+									0, k);
+						com_y += gsl_matrix_get(com,
+									1, k);
 					}
 					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);
+					printf("expanding %i->%i\n", size, size*2);
+		p = itrans_peaksearch_stat_matrix_expand(p, size, size*2);
 						size *= 2;
 					}
-					
+
 				}
 			}
 		}
@@ -452,94 +453,84 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double thresh
 	//printf("ff: ending loop, found %d\n",n);
 
 	*count = n;
-	//printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
-	p = itrans_peaksearch_stat_matrix_expand(p, size, n);
-	//printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
-	
+	if ( n > 0 ) {
+		//printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
+		printf("expandingg %i->%i\n", size, n);
+		p = itrans_peaksearch_stat_matrix_expand(p, size, n);
+		//printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
+	}
+
 	return p;
-	
+
 }
 
 /* Implements the iteration based automatic method
  * returns a gsl_matrix formatted as described in flood-fill */
-static gsl_matrix *itrans_peaksearch_stat_iterate(gsl_matrix *m, unsigned int *count) {
-
+static gsl_matrix *itrans_peaksearch_stat_iterate(gsl_matrix *m,
+							unsigned int *count)
+{
 	int old;
 	int cur;
 	double k;
 	double mean,sigma;
 	gsl_matrix *p;
 	gsl_matrix *kernel;
-	
+
 	old = m->size1*m->size2;
-	
-	//printf("Iterate: starting\n");
-	//printf("Iterate: generating kernel\n");
+
 	kernel = itrans_peaksearch_stat_generate_flat_kernel(1);
-	//printf("Iterate: performing local_sigma_filter\n");
 	itrans_peaksearch_stat_local_sigma_filter(m, 10, 1.);
-	
-	//printf("Iterate: starting loop\n");
+
 	while ( 1 ) {
-	
-		//printf("Iterate: smoothing");
+
 		itrans_peaksearch_stat_apply_kernel(m, kernel);
-		//printf(" (1)");
 		itrans_peaksearch_stat_apply_kernel(m, kernel);
-		//printf(" (2)\n");
 		mean = itrans_peaksearch_stat_image_mean(m);
 		sigma = itrans_peaksearch_stat_image_sigma(m,mean);
-		//printf("Iterate: mean=%lf sigma=%lf\n",mean,sigma);
 		k = (0.5-mean)/sigma;
-		//printf("Iterate: applying sigma_filter(%lf)\n",k);
 		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;
-		
+
 	}
-	
-	gsl_matrix_free(kernel); 
-	//printf("Iterate: finished\n");
+
+	gsl_matrix_free(kernel);
 	*count = cur;
-	
-	return p;
 
+	return p;
 }
 
-ImageFeatureList *itrans_peaksearch_stat(ImageRecord *imagerecord) {
-
-	unsigned int n_reflections;	
+ImageFeatureList *itrans_peaksearch_stat(ImageRecord *imagerecord)
+{
+	unsigned int n_reflections;
 	gsl_matrix *m;
 	gsl_matrix *p;
 	int i;
 	uint16_t *image = imagerecord->image;
 	ImageFeatureList *flist;
-	
+
 	flist = image_feature_list_new();
 
-	m = itrans_peaksearch_stat_createImageMatrix(image, imagerecord->width, imagerecord->height);
+	m = itrans_peaksearch_stat_createImageMatrix(image, imagerecord->width,
+							imagerecord->height);
 	p = itrans_peaksearch_stat_iterate(m, &n_reflections);
-	
+
 	for ( i=0; i<n_reflections; i++ ) {
-	
+
 		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);
-		
+
 	}
-	
+
 	gsl_matrix_free(m);
 	gsl_matrix_free(p);
-	
+
 	return flist;
-	
 }
-