Fix namespace in itrans-stat,{c,h}

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

1 files changed, 82 insertions, 57 deletions

diff --git a/src/itrans-stat.c b/src/itrans-stat.c
index 98988bf..9a876a8 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -1,31 +1,26 @@
 /*
- * peakdetect.c
+ * itrans-stat.c
  *
- * Peakdetection routines, utilities and automation
+ * Peak detection by iterative statistical analysis and processing
  *
  * (c) 2007 Gordon Ball <gfb21@cam.ac.uk>
+ *	    Thomas WHite <taw27@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
 
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
 #include <math.h>
+#include <gsl/gsl_matrix.h>
+
 #include "reflections.h"
 #include "control.h"
 #include "imagedisplay.h"
 #include "utils.h"
-#include "itrans-stat.h"
-
-/*
- * Renormalise a gsl_matrix to 0->1
- * Re-written to use gsl_matrix native functions
- */
-void renormalise(gsl_matrix *m) {
-	double max,min;
-	gsl_matrix_minmax(m,&min,&max);
-	//printf("Renormalising min=%lf max=%lf\n",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
@@ -34,7 +29,7 @@ void renormalise(gsl_matrix *m) {
  * can be done as matrices to save effort
  * Renormalises matrix to 0->1
  */
-gsl_matrix* createImageMatrix(ControlContext *ctx, int16_t *image) {
+static gsl_matrix *itrans_peaksearch_stat_createImageMatrix(ControlContext *ctx, int16_t *image) {
 	gsl_matrix *raw;
 	raw = gsl_matrix_alloc(ctx->width,ctx->height);
 	
@@ -45,7 +40,7 @@ gsl_matrix* createImageMatrix(ControlContext *ctx, int16_t *image) {
 		}
 	}
 	//printf("Created %dx%d matrix\n",ctx->width,ctx->height);
-	renormalise(raw);
+	matrix_renormalise(raw);
 	return raw;
 }
 
@@ -54,7 +49,7 @@ gsl_matrix* createImageMatrix(ControlContext *ctx, int16_t *image) {
 /*
  * Find and return the mean value of the matrix
  */
-double image_mean(gsl_matrix *m) {
+static double itrans_peaksearch_stat_image_mean(gsl_matrix *m) {
 	double mean=0.;
 	int i,j;
 	for (i=0;i<m->size1;i++) {
@@ -69,7 +64,7 @@ double image_mean(gsl_matrix *m) {
  * Return the standard deviation, sigma, of the matrix values
  * \sqrt(\sum((x-mean)^2)/n)
  */
-double image_sigma(gsl_matrix *m, double mean) {
+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++) {
@@ -84,11 +79,11 @@ double image_sigma(gsl_matrix *m, double mean) {
  * Filter all pixels with value < mean + k*sigma to 0
  * Set all matching pixels to 1
  */
-void 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;
-	mean = image_mean(m);
-	sigma = image_sigma(m,mean);
+	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++) {
 			if (gsl_matrix_get(m,i,j) >= mean+k*sigma) {
@@ -105,7 +100,7 @@ void sigma_filter(gsl_matrix *m, double k) {
  *
  * TODO: Use a mask instead of calculating valid points
  */
-double 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,n=0;
 	for (i=x-r;i<=x+r;i++) {
@@ -127,7 +122,7 @@ double circle_mean(gsl_matrix *m, int x, int y, int r, gsl_matrix *mask) {
  *
  * TODO: Use a mask instead of calculating valid points
  */
-double 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,n=0;
 	for (i=x-r;i<=x+r;i++) {
@@ -144,7 +139,7 @@ double circle_sigma(gsl_matrix *m, int x, int y, int r, gsl_matrix *mask, double
 /*
  * Calculate a circular mask to save recalculating it every time
  */
-gsl_matrix* circle_mask(int r) {
+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;
@@ -180,16 +175,16 @@ gsl_matrix* circle_mask(int r) {
  * problem carried over from the OO version where a new object was created by each operation
  */
 
-void 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) {
 	//printf("lsf: starting\n");
 	double mean,sigma;
 	double local;
 	//printf("lsf: generating circle mask\n");
-	gsl_matrix *mask = circle_mask(r);
+	gsl_matrix *mask = itrans_peaksearch_stat_circle_mask(r);
 	gsl_matrix *new;
 	new = gsl_matrix_alloc(m->size1,m->size2);
 	int i,j;
-	int interval = (m->size1-r)/20;
+	//int interval = (m->size1-r)/20;
 	//printf("lsf: starting loop\n");
 	//printf("lsf: ");
 	//for (i=r;i<m->size1-r;i++) {
@@ -198,9 +193,9 @@ void local_sigma_filter(gsl_matrix *m, int r, double k) {
 		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 = circle_mean(m,i,j,r,mask);
+				mean = itrans_peaksearch_stat_circle_mean(m,i,j,r,mask);
 				//printf("lsf: mean=%lf",mean);
-				sigma = 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);
@@ -235,7 +230,7 @@ void local_sigma_filter(gsl_matrix *m, int r, double k) {
  * 
  * Also suffers from self-reference problem
  */
-void apply_kernel(gsl_matrix *m, gsl_matrix *kernel) {
+static void itrans_peaksearch_stat_apply_kernel(gsl_matrix *m, gsl_matrix *kernel) {
 	int size = kernel->size1;
 	int half = (size-1)/2;
 	gsl_matrix *l;
@@ -267,7 +262,7 @@ void apply_kernel(gsl_matrix *m, gsl_matrix *kernel) {
 /*
  * Generate the simplist possible kernel - a flat one
  */
-gsl_matrix* 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)));
@@ -277,7 +272,7 @@ gsl_matrix* generate_flat_kernel(int half) {
 /*
  * expands or contracts a gsl_matrix by copying the columns to a new one
  */
-gsl_matrix* matrix_expand(gsl_matrix *m, int oldsize, int newsize) {
+static gsl_matrix *itrans_peaksearch_stat_matrix_expand(gsl_matrix *m, int oldsize, int newsize) {
 	gsl_matrix *new;
 	
 	//printf("me: %d->%d\n",oldsize,newsize);
@@ -305,7 +300,7 @@ gsl_matrix* matrix_expand(gsl_matrix *m, int oldsize, int newsize) {
  * track of the location of the resized instance
  */
  
-gsl_matrix* 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) {
 	if (i>=0 && i<m->size1) {
 		if (j>=0 && j<m->size2) {
 			if (mask[i+j*m->size1]==0) {
@@ -315,14 +310,14 @@ gsl_matrix* do_ff(int i, int j, int* mask, double threshold, gsl_matrix *m, gsl_
 					gsl_matrix_set(com,1,*com_n,(double)j);
 					*com_n=*com_n+1;
 					if (*com_n == *com_size) {
-						com = 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 = do_ff(i+1,j,mask,threshold,m,com,com_n,com_size);
-					com = do_ff(i-1,j,mask,threshold,m,com,com_n,com_size);
-					com = do_ff(i,j+1,mask,threshold,m,com,com_n,com_size);
-					com = 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);
+					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);
 				}
 			}
 		}
@@ -341,7 +336,7 @@ gsl_matrix* do_ff(int i, int j, int* mask, double threshold, gsl_matrix *m, gsl_
  * Variable count is set to the number of points found
  */
 
-gsl_matrix* floodfill(gsl_matrix *m, double threshold, int *count) {
+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));
@@ -365,7 +360,7 @@ gsl_matrix* floodfill(gsl_matrix *m, double threshold, int *count) {
 					com_x = com_y = 0.;					
 					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 = do_ff(i, j, mask, threshold, m, com, &com_n, &com_size);
+					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);
@@ -378,7 +373,7 @@ gsl_matrix* floodfill(gsl_matrix *m, double threshold, int *count) {
 					gsl_matrix_set(p,1,n,com_y);
 					n++;
 					if (n==size) {
-						p = matrix_expand(p,size,size*2);
+						p = itrans_peaksearch_stat_matrix_expand(p,size,size*2);
 						size *= 2;
 					}
 				}
@@ -388,19 +383,17 @@ gsl_matrix* floodfill(gsl_matrix *m, double threshold, int *count) {
 	//printf("ff: ending loop, found %d\n",n);
 	*count = n;
 	//printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
-	p = 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
- */
-gsl_matrix* iterate(gsl_matrix *m, unsigned int *count) {
-	printf("Iterate: starting\n");
+/* 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 cur;
 	double k;
@@ -408,26 +401,27 @@ gsl_matrix* iterate(gsl_matrix *m, unsigned int *count) {
 	gsl_matrix *p;
 	gsl_matrix *kernel;
 	
+	printf("Iterate: starting\n");
 	//printf("Iterate: generating kernel\n");
-	kernel = generate_flat_kernel(1);
+	kernel = itrans_peaksearch_stat_generate_flat_kernel(1);
 	printf("Iterate: performing local_sigma_filter\n");
-	local_sigma_filter(m,10,1.);
+	itrans_peaksearch_stat_local_sigma_filter(m,10,1.);
 	
 	//printf("Iterate: starting loop\n");
 	while (1) {
 		//printf("Iterate: smoothing");
-		apply_kernel(m,kernel);
+		itrans_peaksearch_stat_apply_kernel(m,kernel);
 		//printf(" (1)");
-		apply_kernel(m,kernel);
+		itrans_peaksearch_stat_apply_kernel(m,kernel);
 		//printf(" (2)\n");
-		mean = image_mean(m);
-		sigma = image_sigma(m,mean);
+		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);
-		sigma_filter(m,k);
+		itrans_peaksearch_stat_sigma_filter(m,k);
 		//printf("Iterate: floodfilling\n");
-		p = floodfill(m,0,&cur);
+		p = itrans_peaksearch_stat_floodfill(m,0,&cur);
 		printf("Iterate: %d points found\n",cur);
 		if (old < 1.05*cur) break;
 		old = cur;
@@ -437,3 +431,34 @@ gsl_matrix* iterate(gsl_matrix *m, unsigned int *count) {
 	*count = cur;
 	return p;
 }
+
+unsigned int itrans_peaksearch_stat(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay) {
+
+	unsigned int n_reflections;	
+	gsl_matrix *m;
+	gsl_matrix *p;
+	int i;
+	double px,py;
+
+	m = itrans_peaksearch_stat_createImageMatrix(ctx, image);
+	p = itrans_peaksearch_stat_iterate(m, &n_reflections);
+	
+	for ( i=0; i<n_reflections; i++ ) {
+	
+		px = gsl_matrix_get(p,0,i);
+		py = gsl_matrix_get(p,1,i);
+		if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
+			reflection_add_from_reciprocal(ctx, (px-ctx->x_centre)*ctx->pixel_size, (py-ctx->y_centre)*ctx->pixel_size,
+								tilt_degrees, 1.0);
+		} else {
+			reflection_add_from_dp(ctx, (px-ctx->x_centre), (py-ctx->y_centre), tilt_degrees, 1.0);
+		}
+		if (ctx->first_image) imagedisplay_mark_point(imagedisplay, (unsigned int)px, (unsigned int)py);
+		
+	}
+	gsl_matrix_free(m);
+	gsl_matrix_free(p);
+	
+	return n_reflections;
+	
+}