From f72af885e7beb97127e4300d635ab156b4336094 Mon Sep 17 00:00:00 2001 From: taw27 Date: Sat, 31 Mar 2007 15:00:56 +0000 Subject: 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 --- src/control.h | 1 + src/displaywindow.c | 1 + src/itrans-stat.c | 330 ++++++++++++++++++++++++++++++---------------------- src/main.c | 4 +- src/mrc.c | 8 +- src/mrc.h | 1 + src/reflections.c | 3 +- src/reflections.h | 7 +- 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 \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 - * Thomas WHite + * Thomas White * * 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; isize1;i++) { - for (j=0; jsize2;j++) { - gsl_matrix_set(raw,i,j,(double)image[i+j*ctx->width]); + for ( i=0; isize1; i++ ) { + for ( j=0; jsize2; 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;isize1;i++) { - for (j=0;jsize2;j++) { + for ( i=0; isize1; i++ ) { + for ( j=0; jsize2; 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;isize1;i++) { - for (j=0;jsize2;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); } } + 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;isize1;i++) { - for (j=0;jsize2;j++) { + + for ( i=0; isize1; i++ ) { + for ( j=0; jsize2; 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;isize1-r;i++) { // for (j=r;jsize2-r;j++) { - for (i=0;isize1;i++) { - for (j=0;jsize2;j++) { - if ((i>=r && isize1-r) && (j>=r && jsize2-r)) { + + 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); + 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;isize1;i++) { - for (j=0;jsize2;j++) { - if ((i>=half && isize1-half) && (j>=half && jsize2-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; 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); l = &lv.matrix; val = 0.; - for (x=0;x%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; jsize1,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 && isize1) { - if (j>=0 && jsize2) { - 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;isize1;i++) { - for (j=0;jsize2;j++) { - if (gsl_matrix_get(m,i,j)>threshold) { - if (mask[i+j*m->size1]==0) { + 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 ) { + //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;ksize1,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 + * * 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; yx_centre = max_x; ctx->y_centre = max_y; + printf("done\n"); for ( i=0; i + * * 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 + * * 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 + * * dtr - Diffraction Tomography Reconstruction * */ +#ifndef REFLECTION_H +#define REFLECTION_H + #ifdef HAVE_CONFIG_H #include #endif -#ifndef REFLECTION_H -#define REFLECTION_H - typedef enum { REFLECTION_NORMAL, /* Measured - expressed as x, y, z position */ REFLECTION_CENTRAL, /* The central beam */ -- cgit v1.2.3