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

7 files changed, 126 insertions, 105 deletions

diff --git a/src/control.h b/src/control.h
index 05a3c5d..2f5f6b9 100644
--- a/src/control.h
+++ b/src/control.h
@@ -34,7 +34,7 @@ typedef enum {
 	PEAKSEARCH_ADAPTIVE_THRESHOLD,
 	PEAKSEARCH_LSQ,
 	PEAKSEARCH_ZAEFFERER,
-	PEAKSEARCH_ITERATIVE
+	PEAKSEARCH_STAT
 } PeakSearchMode;
 
 typedef enum {
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;
+	
+}
diff --git a/src/itrans-stat.h b/src/itrans-stat.h
index a7e0529..ad2f8e3 100644
--- a/src/itrans-stat.h
+++ b/src/itrans-stat.h
@@ -1,19 +1,27 @@
 /*
- * peakdetect.h
+ * itrans-stat.h
+ *
+ * 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
+ *
  */
+ 
+#ifndef ITRANS_STAT_H
+#define ITRANS_STAT_H
 
-#include <gsl/gsl_matrix.h>
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
 
-extern gsl_matrix* createImageMatrix(ControlContext *ctx, int16_t *image);
+#include <stdint.h>
 
-extern void sigma_filter(gsl_matrix *m, double k);
+#include "control.h"
+#include "imagedisplay.h"
 
-extern void local_sigma_filter(gsl_matrix *m, int r, double k);
+unsigned int itrans_peaksearch_stat(int16_t *image, ControlContext *ctx, double tilt_degrees, ImageDisplay *imagedisplay);
 
-extern void apply_kernel(gsl_matrix *m, gsl_matrix *kernel);
-
-extern gsl_matrix* generate_flat_kernel(int half);
-
-extern gsl_matrix* floodfill(gsl_matrix *m, double threshold, int* count);
-
-extern gsl_matrix* iterate(gsl_matrix *m, unsigned int* count);
+#endif /* ITRANS_STAT_H */
diff --git a/src/itrans.c b/src/itrans.c
index cc5941b..d1e43ab 100644
--- a/src/itrans.c
+++ b/src/itrans.c
@@ -25,36 +25,6 @@
 #include "itrans-lsq.h"
 #include "itrans-stat.h"
 
-static unsigned int itrans_peaksearch_iterative(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 = createImageMatrix(ctx, image);
-	p = 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;
-	
-}
-
 void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degrees) {
 	
 	unsigned int n_reflections;
@@ -70,7 +40,7 @@ void itrans_process_image(int16_t *image, ControlContext *ctx, double tilt_degre
 		case PEAKSEARCH_ADAPTIVE_THRESHOLD : n_reflections = itrans_peaksearch_adaptive_threshold(image, ctx, tilt_degrees, imagedisplay); break;
 		case PEAKSEARCH_LSQ : n_reflections = itrans_peaksearch_lsq(image, ctx, tilt_degrees, imagedisplay); break;
 		case PEAKSEARCH_ZAEFFERER : n_reflections = itrans_peaksearch_zaefferer(image, ctx, tilt_degrees, imagedisplay); break;
-		case PEAKSEARCH_ITERATIVE : n_reflections = itrans_peaksearch_iterative(image, ctx, tilt_degrees, imagedisplay); break;
+		case PEAKSEARCH_STAT : n_reflections = itrans_peaksearch_stat(image, ctx, tilt_degrees, imagedisplay); break;
 		default: n_reflections = 0;
 	}
 	
diff --git a/src/main.c b/src/main.c
index 5814e13..dac4379 100644
--- a/src/main.c
+++ b/src/main.c
@@ -49,7 +49,7 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 			case 1 : ctx->psmode = PEAKSEARCH_ADAPTIVE_THRESHOLD; break;
 			case 2 : ctx->psmode = PEAKSEARCH_LSQ; break;
 			case 3 : ctx->psmode = PEAKSEARCH_ZAEFFERER; break;
-			case 4 : ctx->psmode = PEAKSEARCH_ITERATIVE; break;
+			case 4 : ctx->psmode = PEAKSEARCH_STAT; break;
 			default: abort();
 		}
 		
@@ -128,7 +128,7 @@ void main_method_dialog_open(ControlContext *ctx) {
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Adaptive Thresholding");
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Least-Squares Fit");
 	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Zaefferer Gradient Search");
-	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Iterative Statistical");
+	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Iterative Statistical Analysis");
 	gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 3);
 	gtk_table_attach_defaults(GTK_TABLE(table), ctx->combo_peaksearch, 2, 3, 2, 3);
 	
diff --git a/src/utils.c b/src/utils.c
index 1c7ecfb..fdc245f 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -9,6 +9,7 @@
  */
 
 #include <math.h>
+#include <gsl/gsl_matrix.h>
 
 /* Return the MOST POSITIVE of two numbers */
 unsigned int biggest(signed int a, signed int b) {
@@ -54,3 +55,15 @@ double lambda(double V) {
 	return h / sqrt(2*m*e*V*(1+((e*V) / (2*m*c*c))));
 
 }
+
+/* Renormalise a gsl_matrix to 0->1
+ * Re-written to use gsl_matrix native functions */
+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));
+
+}
diff --git a/src/utils.h b/src/utils.h
index b4a0dca..1d5cb3c 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -4,16 +4,20 @@
  * Utility stuff
  *
  * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *	    Gordon Ball <gfb21@cam.ac.uk>
+ *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
  
+#ifndef UTILS_H
+#define UTILS_H
+
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
 
-#ifndef UTILS_H
-#define UTILS_H
+#include <gsl/gsl_matrix.h>
 
 extern unsigned int biggest(signed int a, signed int b);
 extern unsigned int smallest(signed int a, signed int b);
@@ -21,5 +25,6 @@ extern double distance(double x1, double y1, double x2, double y2);
 extern double modulus(double x, double y, double z);
 extern double angle_between(double x1, double y1, double z1, double x2, double y2, double z2);
 extern double lambda(double voltage);
+extern void matrix_renormalise(gsl_matrix *m);
 
 #endif	/* UTILS_H */