Commit ecc4bcba authored by Thomas White's avatar Thomas White
Browse files

indexamajig: Option of getting peaks from the HDF5 file

parent 7ff28769
......@@ -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;
......
......@@ -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 */
......@@ -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;
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment