From 3533748c223433b53318465d2ec16dc72f4d7cdc Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 17 Mar 2009 17:11:24 +0000 Subject: Line wrapping, whitespace, remove debug --- src/itrans-stat.c | 309 ++++++++++++++++++++++++++---------------------------- 1 file 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; isize1; i++ ) { for ( j=0; jsize2; 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;isize1-r;i++) { - // for (j=r;jsize2-r;j++) { - for ( i=0; isize1; i++ ) { for ( j=0; jsize2; 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; isize1; i++ ) { for ( j=0; jsize2; 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%d\n",oldsize,newsize); - + new = gsl_matrix_calloc(3, newsize); - //printf("me: calloc(2,%d)\n",newsize); + for ( j=0; jsize1,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; isize1; i++ ) { for ( j=0; jsize2; 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%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