diff options
author | Thomas White <taw@physics.org> | 2010-10-07 11:03:56 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2012-02-22 15:27:01 +0100 |
commit | ecc4bcba33c84e5f9d8f3909b85c6c81158e99f6 (patch) | |
tree | 703e887eeaa625cb7b9c11727e10da8edee82b3b | |
parent | 7ff287694b7fd5d1fca4f31e552c8d656257c82b (diff) |
indexamajig: Option of getting peaks from the HDF5 file
-rw-r--r-- | src/hdf5-file.c | 78 | ||||
-rw-r--r-- | src/hdf5-file.h | 1 | ||||
-rw-r--r-- | src/indexamajig.c | 50 |
3 files changed, 127 insertions, 2 deletions
diff --git a/src/hdf5-file.c b/src/hdf5-file.c index ed1c5fe9..5022dd44 100644 --- a/src/hdf5-file.c +++ b/src/hdf5-file.c @@ -101,6 +101,84 @@ int hdfile_get_height(struct hdfile *f) } +int get_peaks(struct image *image, struct hdfile *f) +{ + hid_t dh, sh; + hsize_t size[2]; + hsize_t max_size[2]; + int i; + float *buf; + herr_t r; + + dh = H5Dopen(f->fh, "/processing/hitfinder/peakinfo", H5P_DEFAULT); + + if ( dh < 0 ) { + ERROR("No peak list found!\n"); + return 1; + } + + sh = H5Dget_space(dh); + if ( sh < 0 ) { + H5Dclose(dh); + ERROR("Couldn't get dataspace for peak list.\n"); + return 1; + } + + if ( H5Sget_simple_extent_ndims(sh) != 2 ) { + ERROR("Peak list has the wrong dimensionality (%i).\n", + H5Sget_simple_extent_ndims(sh)); + H5Sclose(sh); + H5Dclose(dh); + return 1; + } + + H5Sget_simple_extent_dims(sh, size, max_size); + + if ( size[1] != 3 ) { + H5Sclose(sh); + H5Dclose(dh); + ERROR("Peak list has the wrong dimensions.\n"); + return 1; + } + + buf = malloc(sizeof(float)*size[0]*size[1]); + if ( buf == NULL ) { + H5Sclose(sh); + H5Dclose(dh); + ERROR("Couldn't reserve memory for the peak list.\n"); + return 1; + } + r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf); + if ( r < 0 ) { + ERROR("Couldn't read peak list.\n"); + free(buf); + return 1; + } + + for ( i=0; i<size[0]; i++ ) { + + unsigned int x, y; + float val; + + x = buf[3*i+0]; + y = buf[3*i+1]; + val = buf[3*i+2]; + + STATUS("%i %i %f\n", x, y, val); + + image_add_feature(image->features, x, y, image, val, NULL); + + } + STATUS("Got %i peaks from peak list.\n", (int)size[0]); + + free(buf); + H5Sclose(sh); + H5Dclose(dh); + + return 0; +} + + static void cleanup(hid_t fh) { int n_ids, i; diff --git a/src/hdf5-file.h b/src/hdf5-file.h index 7e662799..145c9cac 100644 --- a/src/hdf5-file.h +++ b/src/hdf5-file.h @@ -41,5 +41,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 int get_peaks(struct image *image, struct hdfile *f); #endif /* HDF5_H */ diff --git a/src/indexamajig.c b/src/indexamajig.c index cebdf9e3..78d5cf57 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -39,6 +39,13 @@ #define MAX_THREADS (96) + +enum { + PEAK_ZAEF, + PEAK_HDF5, +}; + + struct process_args { /* Input */ @@ -67,6 +74,7 @@ struct process_args IndexingPrivate *ipriv; const double *intensities; struct gpu_context *gctx; + int peaks; /* Thread control and output */ pthread_mutex_t control_mutex; /* Protects the scary stuff below */ @@ -100,6 +108,11 @@ static void show_help(const char *s) " -g. --geometry=<file> Get detector geometry from file.\n" " -p, --pdb=<file> PDB file from which to get the unit cell to match.\n" " -x, --prefix=<p> Prefix filenames from input file with <p>.\n" +" --peaks=<method> Use 'method' for finding peaks. Choose from:\n" +" zaef : Use Zaefferer (2000) gradient detection.\n" +" This is the default method.\n" +" hdf5 : Get from /processing/hitfinder/peakinfo\n" +" in the HDF5 file.\n" "\n" "\nWith just the above options, this program does not do much of practical use." "\nYou should also enable some of the following:\n\n" @@ -328,8 +341,19 @@ static void process_image(struct process_args *pargs) memcpy(data_for_measurement, image.data, data_size); } - /* Perform 'fine' peak search */ - search_peaks(&image, pargs->threshold); + switch ( pargs->peaks ) + { + case PEAK_HDF5 : + /* Get peaks from HDF5 */ + if ( get_peaks(&image, hdfile) ) { + ERROR("Failed to get peaks from HDF5 file.\n"); + return; + } + break; + case PEAK_ZAEF : + search_peaks(&image, pargs->threshold); + break; + } /* Get rid of noise-filtered version at this point * - it was strictly for the purposes of peak detection. */ @@ -478,6 +502,8 @@ int main(int argc, char *argv[]) char *intfile = NULL; char *pdb = NULL; char *prefix = NULL; + char *speaks = NULL; + int peaks; int nthreads = 1; pthread_t workers[MAX_THREADS]; struct process_args *worker_args[MAX_THREADS]; @@ -497,6 +523,7 @@ int main(int argc, char *argv[]) {"gpu", 0, &config_gpu, 1}, {"no-index", 0, &config_noindex, 1}, {"dump-peaks", 0, &config_dumpfound, 1}, + {"peaks", 1, NULL, 2}, {"near-bragg", 0, &config_nearbragg, 1}, {"write-drx", 0, &config_writedrx, 1}, {"indexing", 1, NULL, 'z'}, @@ -566,6 +593,10 @@ int main(int argc, char *argv[]) threshold = strtof(optarg, NULL); break; + case 2 : + speaks = strdup(optarg); + break; + case 0 : break; @@ -603,6 +634,20 @@ int main(int argc, char *argv[]) } free(outfile); + if ( speaks == NULL ) { + speaks = strdup("zaef"); + STATUS("You didn't specify a peak detection method.\n"); + STATUS("I'm using 'zaef' for you.\n"); + } + if ( strcmp(speaks, "zaef") == 0 ) { + peaks = PEAK_ZAEF; + } else if ( strcmp(speaks, "hdf5") == 0 ) { + peaks = PEAK_HDF5; + } else { + ERROR("Unrecognised peak detection method '%s'\n", speaks); + return 1; + } + if ( intfile != NULL ) { ReflItemList *items; items = read_reflections(intfile, intensities, @@ -702,6 +747,7 @@ int main(int argc, char *argv[]) worker_args[i] = malloc(sizeof(struct process_args)); worker_args[i]->filename = malloc(1024); worker_args[i]->ofh = ofh; + worker_args[i]->peaks = peaks; worker_active[i] = 0; } |