indexamajig: Option of getting peaks from the HDF5 file

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;
 	}