diff options
Diffstat (limited to 'src/itrans-stat.c')
-rw-r--r-- | src/itrans-stat.c | 139 |
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; + +} |