More fussiness:

	Tidy itrans-stat.c
	Display units for pixel_size
	Display appropriate peak search algorithm for cache files (i.e. none)
	Tweak credits


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

8 files changed, 208 insertions, 147 deletions

diff --git a/src/control.h b/src/control.h
index 2f5f6b9..ed9e128 100644
--- a/src/control.h
+++ b/src/control.h
@@ -30,6 +30,7 @@ typedef enum {
 } FormulationMode;
 
 typedef enum {
+	PEAKSEARCH_NONE,
 	PEAKSEARCH_THRESHOLD,
 	PEAKSEARCH_ADAPTIVE_THRESHOLD,
 	PEAKSEARCH_LSQ,
diff --git a/src/displaywindow.c b/src/displaywindow.c
index f4a8b4a..fa0365a 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -96,6 +96,7 @@ static void displaywindow_about() {
 	gtk_about_dialog_set_comments(GTK_ABOUT_DIALOG(window), "Diffraction Tomography Reconstruction");
 	gtk_about_dialog_set_license(GTK_ABOUT_DIALOG(window), "(c) 2006-2007 Thomas White <taw27@cam.ac.uk>\n"
 								"Virtual trackball (c) Copyright 1993, 1994, Silicon Graphics, Inc.\n"
+								"See Credits for a full list of contributors\n"
 								"\n"
 								"Research funded by:\n"
 								"The Engineering and Physical Sciences Research Council\n"
diff --git a/src/itrans-stat.c b/src/itrans-stat.c
index 9a876a8..c2369ec 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -4,7 +4,7 @@
  * Peak detection by iterative statistical analysis and processing
  *
  * (c) 2007 Gordon Ball <gfb21@cam.ac.uk>
- *	    Thomas WHite <taw27@cam.ac.uk>
+ *	    Thomas White <taw27@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
@@ -22,7 +22,7 @@
 #include "imagedisplay.h"
 #include "utils.h"
 
-/*
+/* 
  * 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
@@ -30,34 +30,36 @@
  * Renormalises matrix to 0->1
  */
 static gsl_matrix *itrans_peaksearch_stat_createImageMatrix(ControlContext *ctx, int16_t *image) {
+
 	gsl_matrix *raw;
+	int i, j;
+
 	raw = gsl_matrix_alloc(ctx->width,ctx->height);
-	
-	int i,j;
-	for ( i=0; i<raw->size1;i++) {
-		for (j=0; j<raw->size2;j++) {
-			gsl_matrix_set(raw,i,j,(double)image[i+j*ctx->width]);
+	for ( i=0; i<raw->size1; i++ ) {
+		for ( j=0; j<raw->size2; j++ ) {
+			gsl_matrix_set(raw, i, j, (double)image[i+j*ctx->width]);
 		}
 	}
-	//printf("Created %dx%d matrix\n",ctx->width,ctx->height);
 	matrix_renormalise(raw);
+	
 	return raw;
+
 }
 
+/* Find and return the mean value of the matrix */
+static double itrans_peaksearch_stat_image_mean(gsl_matrix *m) {
 
+	int i, j;
+	double mean = 0.;
 
-/*
- * Find and return the mean value of the matrix
- */
-static double itrans_peaksearch_stat_image_mean(gsl_matrix *m) {
-	double mean=0.;
-	int i,j;
-	for (i=0;i<m->size1;i++) {
-		for (j=0;j<m->size2;j++) {
+	for ( i=0; i<m->size1; i++ ) {
+		for ( j=0; j<m->size2; j++ ) {
 			mean += gsl_matrix_get(m,i,j);
 		}
 	}
+
 	return mean / (m->size1 * m->size2);
+
 }
 
 /*
@@ -65,34 +67,42 @@ static double itrans_peaksearch_stat_image_mean(gsl_matrix *m) {
  * \sqrt(\sum((x-mean)^2)/n)
  */
 static double itrans_peaksearch_stat_image_sigma(gsl_matrix *m, double mean) {
-	double diff2=0;
+
 	int i,j;
-	for (i=0;i<m->size1;i++) {
-		for (j=0;j<m->size2;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);
 		}
 	}
+	
 	return sqrt(diff2/(m->size1 * m->size2));
+
 }
 
 /*
  * Filter all pixels with value < mean + k*sigma to 0
- * Set all matching pixels to 1
+ * Set all matching pixels to 1 d
  */
 static void itrans_peaksearch_stat_sigma_filter(gsl_matrix *m, double k) {
+
 	double mean,sigma;
 	int i,j;
+
 	mean = itrans_peaksearch_stat_image_mean(m);
 	sigma = itrans_peaksearch_stat_image_sigma(m,mean);
-	for (i=0;i<m->size1;i++) {
-		for (j=0;j<m->size2;j++) {
+
+	for ( i=0; i<m->size1; i++ ) {
+		for ( j=0; j<m->size2; j++ ) {
 			if (gsl_matrix_get(m,i,j) >= mean+k*sigma) {
-				gsl_matrix_set(m,i,j,1.);
+				gsl_matrix_set(m, i, j, 1.);
 			} else {
-				gsl_matrix_set(m,i,j,0.);
+				gsl_matrix_set(m, i, j, 0.);
 			}
 		}
 	}
+	
 }
 
 /*
@@ -101,20 +111,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) {
-	double mean=0.;
-	int i,j,n=0;
-	for (i=x-r;i<=x+r;i++) {
-		for (j=y-r;j<=y+r;j++) {
+
+	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);
+			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;
+	
 }
 
 /*
@@ -123,35 +138,44 @@ 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) {
-	double diff2=0.;
-	int i,j,n=0;
-	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) {
+
+	double diff2 = 0.;
+	int i, j;
+	int n = 0;
+
+	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);
 				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) {
+
 	gsl_matrix *m;
-	m = gsl_matrix_calloc(2*r+1,2*r+1);
 	int i,j;
-	for (i=0;i<2*r+1;i++) {
-		for (j=0;j<2*r+1;j++) {
-			if (sqrt((r-i)*(r-i)+(r-j)*(r-j))<=(double)r) {
-				 gsl_matrix_set(m,i,j,1.);
+
+	m = gsl_matrix_calloc(2*r+1, 2*r+1);
+	for ( i=0; i<2*r+1; i++ ) {
+		for ( j=0; j<2*r+1; j++ ) {
+			if ( sqrt((r-i)*(r-i)+(r-j)*(r-j)) <= r ) {
+				 gsl_matrix_set(m, i, j, 1.);
 			}
 			
 		}
 	}
+	
 	return m;
+	
 }
 
 /*
@@ -174,51 +198,56 @@ static gsl_matrix *itrans_peaksearch_stat_circle_mask(int r) {
  * 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) {
-	//printf("lsf: starting\n");
+
 	double mean,sigma;
 	double local;
-	//printf("lsf: generating circle mask\n");
-	gsl_matrix *mask = itrans_peaksearch_stat_circle_mask(r);
+	gsl_matrix *mask;
 	gsl_matrix *new;
-	new = gsl_matrix_alloc(m->size1,m->size2);
 	int i,j;
-	//int interval = (m->size1-r)/20;
+	//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)) {
+
+	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);
+				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);
+				sigma = itrans_peaksearch_stat_circle_sigma(m, i, j, r, mask, mean);
 				//printf(" sigma=%lf",sigma);
-				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, 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 /= 7.;
 				//printf(" local=%lf\n",local);
-				if (local > mean+k*sigma) {
-					gsl_matrix_set(new,i,j,1.);
+				if ( local > mean+k*sigma ) {
+					gsl_matrix_set(new, i, j, 1.);
 				} else {
-					gsl_matrix_set(new,i,j,0.);
+					gsl_matrix_set(new, i, j, 0.);
 				}
+				
 			} else {
-				gsl_matrix_set(new,i,j,0.);
+				gsl_matrix_set(new, i, j, 0.);
 			}
 		}
 		//if (i % interval == 0) printf(".");
 	}
+	
 	//printf("done\n");
-	gsl_matrix_memcpy(m,new);
+	gsl_matrix_memcpy(m, new);
 	gsl_matrix_free(new);
-}
-
 
+}
 
 /*
  * Apply an arbitary kernel to the image - each point takes the value
@@ -231,67 +260,80 @@ static void itrans_peaksearch_stat_local_sigma_filter(gsl_matrix *m, int r, doub
  * Also suffers from self-reference problem
  */
 static void itrans_peaksearch_stat_apply_kernel(gsl_matrix *m, gsl_matrix *kernel) {
-	int size = kernel->size1;
-	int half = (size-1)/2;
+
+	int size;
+	int half;
 	gsl_matrix *l;
 	gsl_matrix_view lv;
 	gsl_matrix *new;
-	new = gsl_matrix_calloc(m->size1,m->size2);
 	double val;
-	int i,j,x,y;
-	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);
+	int i, j, x, y;
+
+	size = kernel->size1;
+	half = (size-1)/2;
+	new = gsl_matrix_calloc(m->size1, m->size2);
+
+	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);
 				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);
+				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);
 					}
 				}
 				//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
- */
+/* Generate the simplist possible kernel - a flat one */
 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
- */
+/* 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) {
+
 	gsl_matrix *new;
+	int j;
 	
 	//printf("me: %d->%d\n",oldsize,newsize);
 	
-	new = gsl_matrix_calloc(2,newsize);
+	new = gsl_matrix_calloc(2, newsize);
 	//printf("me: calloc(2,%d)\n",newsize);
-	int j;
-	for (j = 0; j < oldsize; j++) {
-		if (j < 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, 0, j, gsl_matrix_get(m, 0, j));
+			gsl_matrix_set(new, 1, j, gsl_matrix_get(m, 1, j));
 		}
 	}
+
 	//printf("me: freeing old matrix\n");
 	gsl_matrix_free(m);
+
 	//printf("me: new s1=%d s2=%d\n",new->size1,new->size2);
-	return new;
 	//printf("me: m s1=%d s2=%d\n",m->size1,m->size2);
+	return new;
+	
 }
 
 /*
@@ -299,30 +341,34 @@ 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) {
-	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) {
+
+	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);
-					gsl_matrix_set(com,0,*com_n,(double)i);
-					gsl_matrix_set(com,1,*com_n,(double)j);
-					*com_n=*com_n+1;
-					if (*com_n == *com_size) {
-						com = itrans_peaksearch_stat_matrix_expand(com,*com_size,*com_size*2);
+					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_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);
+					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);
+					
 				}
 			}
 		}
 	}
+	
 	return com;
+
 }
 
 /*
@@ -335,84 +381,90 @@ static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double
  * coordinates in row 1, which should be of the right length
  * Variable count is set to the number of points found
  */
-
 static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double threshold, int *count) {
-	//printf("ff: starting\n");
+
 	int *mask;
-	mask = calloc(m->size1*m->size2,sizeof(int));
-	
-	int size=32,com_size;
-	int i,j,k,n=0;
+	int size, com_size, i, j, k, n;
 	int com_n;
 	gsl_matrix *p;
 	gsl_matrix *com;
-	p = gsl_matrix_calloc(2,size);
-	
-	double com_x,com_y;
+	double com_x, com_y;
+
+	mask = calloc(m->size1*m->size2, sizeof(int));
+	size = 32;
+	n = 0;
+	p = gsl_matrix_calloc(2, 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) {
+	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 ) {
+				
 					//printf("ff: found starting point (%d,%d)\n",i,j);
-					com_size=32;
-					com_n=0;
+					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(2, 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);
 					//printf("ff: ended floodfill stack\n");
-					for (k=0;k<com_n;k++) {
-						com_x += gsl_matrix_get(com,0,k);
-						com_y += gsl_matrix_get(com,1,k);
+					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 /= com_n;
 					com_y /= com_n;
 					//printf("ff: point CoM (%lf,%lf)\n",com_x,com_y);
-					gsl_matrix_set(p,0,n,com_x);
-					gsl_matrix_set(p,1,n,com_y);
+					gsl_matrix_set(p, 0, n, com_x);
+					gsl_matrix_set(p, 1, n, com_y);
 					n++;
-					if (n==size) {
-						p = itrans_peaksearch_stat_matrix_expand(p,size,size*2);
+					if ( n == size ) {
+						p = itrans_peaksearch_stat_matrix_expand(p, size, size*2);
 						size *= 2;
 					}
+					
 				}
 			}
 		}
 	}
 	//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);
+	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) {
 
-	int old = m->size1*m->size2;
+	int old;
 	int cur;
 	double k;
 	double mean,sigma;
 	gsl_matrix *p;
 	gsl_matrix *kernel;
 	
-	printf("Iterate: starting\n");
+	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: performing local_sigma_filter\n");
+	itrans_peaksearch_stat_local_sigma_filter(m, 10, 1.);
 	
 	//printf("Iterate: starting loop\n");
-	while (1) {
+	while ( 1 ) {
+	
 		//printf("Iterate: smoothing");
-		itrans_peaksearch_stat_apply_kernel(m,kernel);
+		itrans_peaksearch_stat_apply_kernel(m, kernel);
 		//printf(" (1)");
-		itrans_peaksearch_stat_apply_kernel(m,kernel);
+		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);
@@ -421,15 +473,19 @@ static gsl_matrix *itrans_peaksearch_stat_iterate(gsl_matrix *m, unsigned int *c
 		//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);
-		printf("Iterate: %d points found\n",cur);
-		if (old < 1.05*cur) break;
+		p = itrans_peaksearch_stat_floodfill(m, 0, &cur);
+		//printf("Iterate: %d points found\n", cur);
+		if ( old < 1.05*cur ) break;
 		old = cur;
+		
 	}
+	
 	gsl_matrix_free(kernel); 
-	printf("Iterate: finished\n");
+	//printf("Iterate: finished\n");
 	*count = cur;
+	
 	return p;
+
 }
 
 unsigned int itrans_peaksearch_stat(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
diff --git a/src/main.c b/src/main.c
index dac4379..cc302bb 100644
--- a/src/main.c
+++ b/src/main.c
@@ -50,7 +50,7 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 			case 2 : ctx->psmode = PEAKSEARCH_LSQ; break;
 			case 3 : ctx->psmode = PEAKSEARCH_ZAEFFERER; break;
 			case 4 : ctx->psmode = PEAKSEARCH_STAT; break;
-			default: abort();
+			default: ctx->psmode = PEAKSEARCH_NONE; break;	/* This happens when reading from a cache file */
 		}
 		
 		gtk_widget_destroy(method_window);
@@ -133,7 +133,9 @@ void main_method_dialog_open(ControlContext *ctx) {
 	gtk_table_attach_defaults(GTK_TABLE(table), ctx->combo_peaksearch, 2, 3, 2, 3);
 	
 	if ( ctx->inputfiletype == INPUT_CACHE ) {
+		gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Get from cache file");
 		gtk_widget_set_sensitive(GTK_WIDGET(ctx->combo_peaksearch), FALSE);
+		gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 5);
 	}
 
 	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(table), TRUE, TRUE, 5);
diff --git a/src/mrc.c b/src/mrc.c
index a01809c..f09c79e 100644
--- a/src/mrc.c
+++ b/src/mrc.c
@@ -4,6 +4,7 @@
  * Read the MRC tomography format
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
@@ -55,7 +56,6 @@ int mrc_read(ControlContext *ctx) {
 	
 	/* Read all extended headers, one by one */
 	extsize = 4*mrc.numintegers + 4*mrc.numfloats;
-	printf("extsize=%d numintegers=%d numfloats=%d\n",extsize,mrc.numintegers,mrc.numfloats);
 	if ( extsize > sizeof(MRCExtHeader) ) {
 		fclose(fh);
 		fprintf(stderr, "MR: MRC extended header is too big - tweak mrc.h\n");
@@ -66,14 +66,14 @@ int mrc_read(ControlContext *ctx) {
 	}
 	
 	pixel_size = ext[0].pixel_size;
-	printf("pixel_size=%f\n", pixel_size);
+	printf("pixel_size = %f m^-1\n", pixel_size);
 	ctx->reflectionctx = reflection_init();
 	ctx->fmode = FORMULATION_PIXELSIZE;
 	ctx->first_image = 1;
 	ctx->width = mrc.nx;
 	ctx->height = mrc.ny;
 	
-	printf("Judging centre...\n");
+	printf("Judging centre...");
 	image_total = malloc(mrc.ny * mrc.nx * sizeof(int16_t));
 	for ( y=0; y<mrc.ny; y++ ) {
 		for ( x=0; x<mrc.nx; x++ ) {
@@ -102,9 +102,9 @@ int mrc_read(ControlContext *ctx) {
 		}
 	}
 	imagedisplay_mark_point(sum_id, max_x, max_y);
-	printf("max_x=%d max_y=%d\n",max_x,max_y);
 	ctx->x_centre = max_x;
 	ctx->y_centre = max_y;
+	printf("done\n");
 	
 	for ( i=0; i<mrc.nz; i++ ) {
 	
diff --git a/src/mrc.h b/src/mrc.h
index 8633f05..bdef3f8 100644
--- a/src/mrc.h
+++ b/src/mrc.h
@@ -4,6 +4,7 @@
  * Read the MRC tomography format
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
diff --git a/src/reflections.c b/src/reflections.c
index 65e88a6..4ecc268 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -4,6 +4,7 @@
  * Data structures in 3D space
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
@@ -30,8 +31,6 @@ static void reflection_addfirst(ReflectionContext *reflectionctx) {
 	reflectionctx->last_reflection = reflectionctx->reflections;
 }
 
-/* Size of the volume (arbitrary units - m-1 in DTR), number of reflections across.  The centre is defined to be at nx/2, ny/2, nz/2.
-	This defines the three-dimensional grid of reflections */
 ReflectionContext *reflection_init() {
 
 	ReflectionContext *reflectionctx = malloc(sizeof(ReflectionContext));
diff --git a/src/reflections.h b/src/reflections.h
index ac08133..624c1ce 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -4,17 +4,18 @@
  * Data structures in 3D space
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
  
+#ifndef REFLECTION_H
+#define REFLECTION_H
+
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
 
-#ifndef REFLECTION_H
-#define REFLECTION_H
-
 typedef enum {
 	REFLECTION_NORMAL,	/* Measured - expressed as x, y, z position */
 	REFLECTION_CENTRAL,	/* The central beam */