From bd26d5745269594647ec79f64fdfb8e750891672 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 8 Jan 2010 15:53:17 +0100 Subject: Zaefferer gradient search --- src/Makefile.am | 2 +- src/dirax.c | 99 +++++++++++++++++++++++++++++++++++++++++++++++------ src/displaywindow.c | 35 ++++++++++++++++++- src/hdf5-file.c | 22 +++++++++--- src/hdf5-file.h | 1 + src/indexamajig.c | 4 +-- src/render.c | 27 +++++++++++++++ src/utils.h | 4 +++ 8 files changed, 174 insertions(+), 20 deletions(-) diff --git a/src/Makefile.am b/src/Makefile.am index 64567c2c..839f7450 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -20,7 +20,7 @@ indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c dirax.c cell.c indexamajig_LDADD = @LIBS@ if HAVE_GTK -hdfsee_SOURCES = hdfsee.c displaywindow.c render.c hdf5-file.c utils.c +hdfsee_SOURCES = hdfsee.c displaywindow.c render.c hdf5-file.c utils.c image.c hdfsee_LDADD = @LIBS@ endif diff --git a/src/dirax.c b/src/dirax.c index c9c5f010..c1794764 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -32,6 +32,7 @@ #include "image.h" #include "dirax.h" +#include "utils.h" typedef enum { @@ -371,10 +372,13 @@ static int map_position(struct image *image, double x, double y, } +#define PEAK_WINDOW_SIZE (10) + static void search_peaks(struct image *image) { FILE *fh; - int x, y; + int x, y, width; + int16_t *data; fh = fopen("xfel.drx", "w"); if ( !fh ) { @@ -383,23 +387,96 @@ static void search_peaks(struct image *image) } fprintf(fh, "%f\n", 0.5); /* Lie about the wavelength. */ - for ( x=0; xwidth; x++ ) { - for ( y=0; yheight; y++ ) { + data = image->data; + width = image->width; + + for ( x=1; xwidth-1; x++ ) { + for ( y=1; yheight-1; y++ ) { + + double dx1, dx2, dy1, dy2; + double dxs, dys; + double grad; + int mask_x, mask_y; + int sx, sy; + double max; + unsigned int did_something = 1; + + /* Get gradients */ + dx1 = data[x+width*y] - data[(x+1)+width*y]; + dx2 = data[(x-1)+width*y] - data[x+width*y]; + dy1 = data[x+width*y] - data[(x+1)+width*(y+1)]; + dy2 = data[x+width*(y-1)] - data[x+width*y]; + + /* Average gradient measurements from both sides */ + dxs = ((dx1*dx1) + (dx2*dx2)) / 2; + dys = ((dy1*dy1) + (dy2*dy2)) / 2; + + /* Calculate overall gradient */ + grad = dxs + dys; - int val; + if ( grad < 2000000 ) continue; - val = image->data[y+image->height*(image->width-1-x)]; + mask_x = x; + mask_y = y; - if ( val > 1000 ) { + while ( (did_something) + && (distance(mask_x, mask_y, x, y)<50) ) { - double rx, ry, rz; + max = data[mask_x+width*mask_y]; + did_something = 0; + + for ( sy=biggest(mask_y-PEAK_WINDOW_SIZE/2, 0); + syheight); + sy++ ) { + for ( sx=biggest(mask_x-PEAK_WINDOW_SIZE/2, 0); + sx max ) { + max = data[sx+width*sy]; + mask_x = sx; + mask_y = sy; + did_something = 1; + } + + } + } - /* Map and record reflection */ - map_position(image, x, y, &rx, &ry, &rz); - fprintf(fh, "%10f %10f %10f %8f\n", - rx/1e10, ry/1e10, rz/1e10, 1.0); } + if ( !did_something ) { + + //double d; + //int idx; + + assert(mask_xwidth); + assert(mask_yheight); + assert(mask_x>=0); + assert(mask_y>=0); + + if ( distance(mask_x, mask_y, x, y) > 50.0 ) continue; + + /* Check for a feature at exactly the + * same coordinates */ + //image_feature_closest(flist, mask_x, mask_y, + // &d, &idx); + + //if ( d > 1.0 ) { + + double rx = 0.0; + double ry = 0.0; + double rz = 0.0; + + /* Map and record reflection */ + printf("%i %i\n", x, y); + map_position(image, x, y, &rx, &ry, &rz); + fprintf(fh, "%10f %10f %10f %8f\n", + rx/1e10, ry/1e10, rz/1e10, 1.0); + //} + + } + + } } diff --git a/src/displaywindow.c b/src/displaywindow.c index 84f89347..9a65fb16 100644 --- a/src/displaywindow.c +++ b/src/displaywindow.c @@ -26,6 +26,7 @@ #include "render.h" #include "hdf5-file.h" #include "hdfsee.h" +#include "utils.h" #define INITIAL_BINNING 2 @@ -399,6 +400,37 @@ static gint displaywindow_set_boostint(GtkWidget *widget, DisplayWindow *dw) } +static void load_features_from_file(struct image *image, const char *filename) +{ + FILE *fh; + char *rval; + + fh = fopen(filename, "r"); + if ( fh == NULL ) return; + + if ( image->features != NULL ) { + image_feature_list_free(image->features); + } + image->features = image_feature_list_new(); + + do { + char line[1024]; + int x, y; + int r; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + chomp(line); + + r = sscanf(line, "%i %i", &x, &y); + if ( r != 2 ) continue; + + image_add_feature(image->features, x, y, image, 1.0); + + } while ( rval != NULL ); +} + + static gint displaywindow_peaklist_response(GtkWidget *d, gint response, DisplayWindow *dw) { @@ -408,7 +440,8 @@ static gint displaywindow_peaklist_response(GtkWidget *d, gint response, filename = gtk_file_chooser_get_filename(GTK_FILE_CHOOSER(d)); - load_features_from_file(dw->image); + load_features_from_file(hdfile_get_image(dw->hdfile), filename); + displaywindow_update(dw); g_free(filename); diff --git a/src/hdf5-file.c b/src/hdf5-file.c index 0e9453cf..509147a0 100644 --- a/src/hdf5-file.c +++ b/src/hdf5-file.c @@ -45,6 +45,7 @@ struct hdfile *hdfile_open(const char *filename) f = malloc(sizeof(struct hdfile)); if ( f == NULL ) return NULL; + f->image = NULL; f->fh = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT); if ( f->fh < 0 ) { @@ -149,18 +150,29 @@ int16_t *hdfile_get_image_binned(struct hdfile *f, int binning, int16_t *max) struct image *image; int16_t *data; - image = malloc(sizeof(struct image)); - if ( image == NULL ) return NULL; + if ( f->image == NULL ) { - hdf5_read(f, image); - f->image = image; + image = malloc(sizeof(struct image)); + if ( image == NULL ) return NULL; + image->features = NULL; - data = hdfile_bin(image->data, f->nx, f->ny, binning, max); + hdf5_read(f, image); + f->image = image; + + } + + data = hdfile_bin(f->image->data, f->nx, f->ny, binning, max); return data; } +struct image *hdfile_get_image(struct hdfile *f) +{ + return f->image; +} + + int hdfile_get_unbinned_value(struct hdfile *f, int x, int y, int16_t *val) { if ( (x>=f->image->width) || (y>=f->image->height) ) { diff --git a/src/hdf5-file.h b/src/hdf5-file.h index a212f677..f0954fbd 100644 --- a/src/hdf5-file.h +++ b/src/hdf5-file.h @@ -42,5 +42,6 @@ extern int hdfile_set_first_image(struct hdfile *f, const char *group); extern void hdfile_close(struct hdfile *f); extern char *hdfile_get_string_value(struct hdfile *f, const char *name); +extern struct image *hdfile_get_image(struct hdfile *f); #endif /* HDF5_H */ diff --git a/src/indexamajig.c b/src/indexamajig.c index 952433dd..6f09df30 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -176,7 +176,7 @@ static int image_fom(struct image *image) n_valid = 0; for ( i=0; i 0 ) { - STATUS("Hit: %s\n", line); + printf("%s\n", line); index_pattern(&image); diff --git a/src/render.c b/src/render.c index 4de441f3..ce44610c 100644 --- a/src/render.c +++ b/src/render.c @@ -89,6 +89,31 @@ static void render_free_data(guchar *data, gpointer p) free(data); } + +static void show_marked_features(struct image *image, guchar *data, + int w, int h, int binning) +{ + int i; + + if ( image->features == NULL ) return; + + for ( i=0; ifeatures); i++ ) { + + struct imagefeature *f; + int x, y; + + f = image_get_feature(image->features, i); + + x = f->x; y = f->y; + + x /= binning; + y /= binning; + + data[3*( x+w*(h-1-y) )+0] = 255; + } +} + + /* Return a pixbuf containing a rendered version of the image after binning. * This pixbuf might be scaled later - hopefully mostly in a downward * direction. */ @@ -144,6 +169,8 @@ GdkPixbuf *render_get_image(struct hdfile *hdfile, int binning, int boostint, } } + show_marked_features(hdfile_get_image(hdfile), data, w, h, binning); + /* Finished with this */ free(hdr); diff --git a/src/utils.h b/src/utils.h index 4f1dc86c..9aa07c8c 100644 --- a/src/utils.h +++ b/src/utils.h @@ -110,6 +110,10 @@ static inline double distance3d(double x1, double y1, double z1, #define is_odd(a) ((a)%2==1) +#define biggest(a,b) ((a>b) ? (a) : (b)) +#define smallest(a,b) ((a