diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-11-30 17:53:27 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-11-30 17:53:27 +0000 |
commit | cf7a9b848cfa0ec4b649f54100b6ff92e5553232 (patch) | |
tree | 6007c5f513fc8fa57a6c78690c3ed1bde927bb6c | |
parent | 627e71291eb879b77a97e1c429cfb8c6181f7179 (diff) |
Attempt to make itrans-stat quantitative
Size of ImageDisplayMark circles depend on intensity of reflection
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@216 bf6ca9ba-c028-0410-8290-897cf20841d1
-rw-r--r-- | src/imagedisplay.c | 24 | ||||
-rw-r--r-- | src/imagedisplay.h | 3 | ||||
-rw-r--r-- | src/itrans-stat.c | 41 | ||||
-rw-r--r-- | src/refine.c | 8 | ||||
-rw-r--r-- | src/reproject.c | 6 |
5 files changed, 57 insertions, 25 deletions
diff --git a/src/imagedisplay.c b/src/imagedisplay.c index 1caaf9e..9c3f6cd 100644 --- a/src/imagedisplay.c +++ b/src/imagedisplay.c @@ -173,6 +173,7 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even double scale, xoffs, yoffs; ImageDisplayMark *cur; + double max; xoffs = ((double)imagedisplay->drawingarea_width - imagedisplay->view_width) / 2; yoffs = ((double)imagedisplay->drawingarea_height - imagedisplay->view_height) / 2; @@ -211,6 +212,15 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even imagedisplay->imagerecord.x_centre * scale, imagedisplay->imagerecord.y_centre * scale + 10); } + + /* Find the maximum intensity */ + cur = imagedisplay->marks; + max = 0.0; + while ( cur ) { + if ( cur->weight > max ) max = cur->weight; + if ( cur->weight < 0.0 ) printf("ID: Warning: ImageDisplayMark with negative weight\n"); + cur = cur->next; + } cur = imagedisplay->marks; while ( cur ) { @@ -229,10 +239,14 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even if ( (cur->type == IMAGEDISPLAY_MARK_CIRCLE_1) || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_2) || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_3) ) { + double r; + + r = 20 * (cur->weight/max); + gdk_draw_arc(drawingarea->window, gc, FALSE, - xoffs + cur->x*scale - 5, - yoffs + imagedisplay->view_height-1-cur->y*scale - 5, - 11, 11, 0, 64*360); + xoffs + cur->x*scale - r, + yoffs + imagedisplay->view_height-1-cur->y*scale - r, + 2*r, 2*r, 0, 64*360); } else if ( (cur->type == IMAGEDISPLAY_MARK_LINE_1) || (cur->type == IMAGEDISPLAY_MARK_LINE_2) ) { @@ -344,13 +358,14 @@ ImageDisplay *imagedisplay_open(ImageRecord image, const char *title, ImageDispl return imagedisplay_open_with_message(image, title, NULL, flags, NULL, NULL); } -void imagedisplay_add_mark(ImageDisplay *imagedisplay, double x, double y, ImageDisplayMarkType type) { +void imagedisplay_add_mark(ImageDisplay *imagedisplay, double x, double y, ImageDisplayMarkType type, double weight) { ImageDisplayMark *new; new = malloc(sizeof(ImageDisplayMark)); new->x = x; new->y = y; new->type = type; + new->weight = weight; new->next = NULL; if ( !imagedisplay->marks ) { @@ -373,6 +388,7 @@ void imagedisplay_add_line(ImageDisplay *imagedisplay, double x1, double y1, dou new->x = x1; new->y = y1; new->x2 = x2; new->y2 = y2; new->type = type; + new->weight = 1.0; /* This field makes little sense for a line */ new->next = NULL; if ( !imagedisplay->marks ) { diff --git a/src/imagedisplay.h b/src/imagedisplay.h index 7b16dc2..681de79 100644 --- a/src/imagedisplay.h +++ b/src/imagedisplay.h @@ -44,6 +44,7 @@ typedef struct struct_imagedisplaymark { double x2; double y2; ImageDisplayMarkType type; + double weight; /* The 'intensity' if this is a peak */ struct struct_imagedisplaymark *next; @@ -81,7 +82,7 @@ typedef struct imagedisplay_struct { extern ImageDisplay *imagedisplay_open(ImageRecord image, const char *title, ImageDisplayFlags flags); extern ImageDisplay *imagedisplay_open_with_message(ImageRecord image, const char *title, const char *message, ImageDisplayFlags flags, GCallback mouse_click_func, gpointer callback_data); -extern void imagedisplay_add_mark(ImageDisplay *imagedisplay, double x, double y, ImageDisplayMarkType type); +extern void imagedisplay_add_mark(ImageDisplay *imagedisplay, double x, double y, ImageDisplayMarkType type, double weight); extern void imagedisplay_add_line(ImageDisplay *imagedisplay, double x1, double y1, double x2, double y2, ImageDisplayMarkType type); extern void imagedisplay_force_redraw(ImageDisplay *imagedisplay); extern void imagedisplay_put_data(ImageDisplay *imagedisplay, ImageRecord imagerecord); diff --git a/src/itrans-stat.c b/src/itrans-stat.c index 67ed96b..fe08a3b 100644 --- a/src/itrans-stat.c +++ b/src/itrans-stat.c @@ -329,13 +329,14 @@ static gsl_matrix *itrans_peaksearch_stat_matrix_expand(gsl_matrix *m, int oldsi //printf("me: %d->%d\n",oldsize,newsize); - new = gsl_matrix_calloc(2, newsize); + new = gsl_matrix_calloc(3, newsize); //printf("me: calloc(2,%d)\n",newsize); for ( j=0; j<oldsize; j++) { if ( j < newsize ) { //printf("me: copying col %d\n",j); gsl_matrix_set(new, 0, j, gsl_matrix_get(m, 0, j)); gsl_matrix_set(new, 1, j, gsl_matrix_get(m, 1, j)); + gsl_matrix_set(new, 2, j, gsl_matrix_get(m, 2, j)); } } @@ -353,7 +354,7 @@ 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) { +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) ) { @@ -361,6 +362,7 @@ static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double 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; @@ -369,10 +371,10 @@ static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double *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); + 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); } } @@ -390,7 +392,8 @@ static gsl_matrix *itrans_peaksearch_stat_do_ff(int i, int j, int* mask, double * dependent on implementation (default 4MiB stack for unix?) * Implements a crude variable-sized array method which hopefully works * Returns a gsl_matrix with x coordinates in row 0 and y - * coordinates in row 1, which should be of the right length + * coordinates in row 1, which should be of the right length. + * 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) { @@ -405,7 +408,7 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double thresh mask = calloc(m->size1*m->size2, sizeof(int)); size = 32; n = 0; - p = gsl_matrix_calloc(2, size); + p = gsl_matrix_calloc(3, size); //printf("ff: starting loop\n"); for ( i=0; i<m->size1; i++ ) { @@ -413,13 +416,16 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double thresh 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(2, com_size); //this is going to hold the points found for this location + 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 = itrans_peaksearch_stat_do_ff(i, j, mask, threshold, m, com, &com_n, &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"); for ( k=0; k<com_n; k++ ) { com_x += gsl_matrix_get(com, 0, k); @@ -428,8 +434,11 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, double thresh com_x /= com_n; com_y /= com_n; //printf("ff: point CoM (%lf,%lf)\n",com_x,com_y); + + /* Now add it to the list */ gsl_matrix_set(p, 0, n, com_x); gsl_matrix_set(p, 1, n, com_y); + gsl_matrix_set(p, 2, n, val); n++; if ( n == size ) { p = itrans_peaksearch_stat_matrix_expand(p, size, size*2); @@ -486,6 +495,8 @@ static gsl_matrix *itrans_peaksearch_stat_iterate(gsl_matrix *m, unsigned int *c 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; @@ -506,7 +517,6 @@ ImageFeatureList *itrans_peaksearch_stat(ImageRecord *imagerecord) { gsl_matrix *m; gsl_matrix *p; int i; - double px,py; uint16_t *image = imagerecord->image; ImageFeatureList *flist; @@ -517,9 +527,12 @@ ImageFeatureList *itrans_peaksearch_stat(ImageRecord *imagerecord) { for ( i=0; i<n_reflections; i++ ) { - px = gsl_matrix_get(p,0,i); - py = gsl_matrix_get(p,1,i); - image_add_feature(flist, px, py, imagerecord, 1.0); + double px, py, pi; + + px = gsl_matrix_get(p, 0, i); + py = gsl_matrix_get(p, 1, i); + pi = gsl_matrix_get(p, 2, i); + image_add_feature(flist, px, py, imagerecord, pi); } diff --git a/src/refine.c b/src/refine.c index 6b6eefc..8b6c716 100644 --- a/src/refine.c +++ b/src/refine.c @@ -172,17 +172,17 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio /* Update the cell */ mapping_rotate(dax, day, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); - printf("a changed by %f %f %f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, + printf("a changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, 100*dlx/cell->a.x, 100*dly/cell->a.y, 100*dlz/cell->a.z); cell->a.x += dlx; cell->a.y += dly; cell->a.z += dlz; mapping_rotate(dbx, dby, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); - printf("b changed by %f %f %f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, + printf("b changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, 100*dlx/cell->b.x, 100*dly/cell->b.y, 100*dlz/cell->b.z); cell->b.x += dlx; cell->b.y += dly; cell->b.z += dlz; mapping_rotate(dcx, dcy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); - printf("c changed by %f %f %f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, + printf("c changed by %+7.5f %+7.5f %+7.5f nm^-1 in reciprocal space (%+10.5f %+10.5f %+10.5f %%)\n", dlx/1e9, dly/1e9, dlz/1e9, 100*dlx/cell->c.x, 100*dly/cell->c.y, 100*dlz/cell->c.z); cell->c.x += dlx; cell->c.y += dly; cell->c.z += dlz; @@ -334,7 +334,7 @@ static gint refine_dialog_fit_image(GtkWidget *step_button, ControlContext *ctx) displaywindow_update(ctx->dw); if ( fitted ) { - imagedisplay_add_mark(ctx->reproject_id, fitted->x,fitted->y, IMAGEDISPLAY_MARK_CIRCLE_3); + imagedisplay_add_mark(ctx->reproject_id, fitted->x,fitted->y, IMAGEDISPLAY_MARK_CIRCLE_3, fitted->intensity); } } else { diff --git a/src/reproject.c b/src/reproject.c index 1fd061a..2bb33dc 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -217,13 +217,15 @@ static void reproject_mark_peaks(ControlContext *ctx) { image->rflist = reproject_get_reflections(image, ctx->cell_lattice); } for ( j=0; j<image->rflist->n_features; j++ ) { - imagedisplay_add_mark(ctx->reproject_id, image->rflist->features[j].x, image->rflist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_1); + imagedisplay_add_mark(ctx->reproject_id, image->rflist->features[j].x, image->rflist->features[j].y, + IMAGEDISPLAY_MARK_CIRCLE_1, image->rflist->features[j].intensity); } /* Now draw the original measured peaks */ flist = image->features; for ( j=0; j<flist->n_features; j++ ) { - imagedisplay_add_mark(ctx->reproject_id, flist->features[j].x, flist->features[j].y, IMAGEDISPLAY_MARK_CIRCLE_2); + imagedisplay_add_mark(ctx->reproject_id, flist->features[j].x, flist->features[j].y, + IMAGEDISPLAY_MARK_CIRCLE_2, image->rflist->features[j].intensity); } /* Now connect partners */ |