indexamajig.c 28.1 KB
Newer Older
Thomas White's avatar
Thomas White committed
1
/*
2
 * indexamajig.c
Thomas White's avatar
Thomas White committed
3
 *
Thomas White's avatar
Thomas White committed
4
 * Index patterns, output hkl+intensity etc.
Thomas White's avatar
Thomas White committed
5
 *
Thomas White's avatar
Thomas White committed
6
 * Copyright © 2012 Thomas White <taw@physics.org>
Thomas White's avatar
Thomas White committed
7
 *
Thomas White's avatar
Thomas White committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 * This file is part of CrystFEL.
 *
 * CrystFEL is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CrystFEL is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
Thomas White's avatar
Thomas White committed
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 *
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <getopt.h>
Thomas White's avatar
Thomas White committed
36
#include <hdf5.h>
Thomas White's avatar
Thomas White committed
37
#include <gsl/gsl_errno.h>
38
#include <pthread.h>
39
40

#ifdef HAVE_CLOCK_GETTIME
41
#include <time.h>
42
43
44
#else
#include <sys/time.h>
#endif
Thomas White's avatar
Thomas White committed
45
46
47

#include "utils.h"
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
48
#include "index.h"
49
#include "peaks.h"
50
#include "detector.h"
Thomas White's avatar
Thomas White committed
51
#include "filters.h"
52
#include "thread-pool.h"
53
#include "beam-parameters.h"
54
#include "geometry.h"
Thomas White's avatar
Thomas White committed
55
#include "stream.h"
Thomas White's avatar
Thomas White committed
56
#include "reflist-utils.h"
57

58

Thomas White's avatar
Thomas White committed
59
/* Write statistics at APPROXIMATELY this interval */
60
61
62
#define STATS_EVERY_N_SECONDS (5)


63
64
65
66
67
68
enum {
	PEAK_ZAEF,
	PEAK_HDF5,
};


69
70
/* Information about the indexing process which is common to all patterns */
struct static_index_args
71
72
73
74
75
{
	UnitCell *cell;
	int config_cmfilter;
	int config_noisefilter;
	int config_verbose;
Thomas White's avatar
Thomas White committed
76
	int stream_flags;         /* What goes into the output? */
Thomas White's avatar
Thomas White committed
77
	int config_polar;
78
	int config_satcorr;
79
	int config_closer;
80
	int config_insane;
81
	int config_bgsub;
82
	float threshold;
83
	float min_gradient;
84
	float min_snr;
85
	double min_int_snr;
Thomas White's avatar
Thomas White committed
86
	struct detector *det;
Thomas White's avatar
Thomas White committed
87
88
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
Thomas White's avatar
Thomas White committed
89
	int peaks;                /* Peak detection method */
90
	int cellr;
91
	float tols[4];
92
	struct beam_params *beam;
Thomas White's avatar
Thomas White committed
93
	const char *element;
94
	const char *hdf5_peak_path;
95

96
97
98
	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
99
	const struct copy_hdf5_field *copyme;
100
101
102
};


103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
/* Information about the indexing process for one pattern */
struct index_args
{
	/* "Input" */
	char *filename;
	struct static_index_args static_args;

	/* "Output" */
	int indexable;
};


/* Information needed to choose the next task and dispatch it */
struct queue_args
{
	FILE *fh;
	char *prefix;
120
	int config_basename;
121
122
	struct static_index_args static_args;

123
	char *use_this_one_instead;
124
125
126
127
128
129

	int n_indexable;
	int n_processed;
	int n_indexable_last_stats;
	int n_processed_last_stats;
	int t_last_stats;
130
131
132
};


Thomas White's avatar
Thomas White committed
133
134
135
136
137
138
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
"Process and index FEL diffraction images.\n"
"\n"
139
" -h, --help               Display this help message.\n"
Thomas White's avatar
Thomas White committed
140
"\n"
141
" -i, --input=<filename>   Specify file containing list of images to process.\n"
Thomas White's avatar
Thomas White committed
142
"                           '-' means stdin, which is the default.\n"
Thomas White's avatar
Thomas White committed
143
144
" -o, --output=<filename>  Write output stream to this file. '-' for stdout.\n"
"                           Default: indexamajig.stream\n"
Thomas White's avatar
Thomas White committed
145
"\n"
Thomas White's avatar
Thomas White committed
146
147
148
149
150
"     --indexing=<methods> Use 'methods' for indexing.  Provide one or more\n"
"                           methods separated by commas.  Choose from:\n"
"                            none     : no indexing (default)\n"
"                            dirax    : invoke DirAx\n"
"                            mosflm   : invoke MOSFLM (DPS)\n"
Thomas White's avatar
Thomas White committed
151
"                            reax     : DPS algorithm with known unit cell\n"
152
" -g. --geometry=<file>    Get detector geometry from file.\n"
153
154
155
" -b, --beam=<file>        Get beam parameters from file (provides nominal\n"
"                           wavelength value if no per-shot value is found in\n"
"                           the HDF5 files.\n"
156
" -p, --pdb=<file>         PDB file from which to get the unit cell to match.\n"
Thomas White's avatar
Thomas White committed
157
"                           Default: 'molecule.pdb'.\n"
158
"     --basename           Remove the directory parts of the filenames.\n"
159
" -x, --prefix=<p>         Prefix filenames from input file with <p>.\n"
160
161
162
"     --peaks=<method>     Use 'method' for finding peaks.  Choose from:\n"
"                           zaef  : Use Zaefferer (2000) gradient detection.\n"
"                                    This is the default method.\n"
163
164
165
"                           hdf5  : Get from a table in HDF5 file.\n"
"     --hdf5-peaks=<p>     Find peaks table in HDF5 file here.\n"
"                           Default: /processing/hitfinder/peakinfo\n"
Thomas White's avatar
Thomas White committed
166
167
"\n\n"
"You can control what information is included in the output stream using\n"
Thomas White's avatar
Thomas White committed
168
"' --record=<flag1>,<flag2>,<flag3>' and so on.  Possible flags are:\n\n"
169
170
" integrated        Include a list of reflection intensities, produced by\n"
"                    integrating around predicted peak locations.\n"
Thomas White's avatar
Thomas White committed
171
"\n"
172
173
" peaks             Include peak locations and intensities from the peak\n"
"                    search.\n"
Thomas White's avatar
Thomas White committed
174
"\n"
175
" peaksifindexed    As 'peaks', but only if the pattern could be indexed.\n"
Thomas White's avatar
Thomas White committed
176
"\n"
177
178
" peaksifnotindexed As 'peaks', but only if the pattern could NOT be indexed.\n"
"\n\n"
Thomas White's avatar
Thomas White committed
179
"The default is '--record=integrated'.\n"
Thomas White's avatar
Thomas White committed
180
181
"\n\n"
"For more control over the process, you might need:\n\n"
Thomas White's avatar
Thomas White committed
182
183
184
185
186
"     --cell-reduction=<m> Use <m> as the cell reduction method. Choose from:\n"
"                           none    : no matching, just use the raw cell.\n"
"                           reduce  : full cell reduction.\n"
"                           compare : match by at most changing the order of\n"
"                                     the indices.\n"
187
"                           compare_ab : compare 'a' and 'b' lengths only.\n"
188
189
190
"     --tolerance=<a,b,c,angl>  Set the tolerance for a,b,c axis (in %%)\n"
"                                and for the angles (in deg) when reducing\n"
"                                or comparing (default is 5%% and 1.5deg)\n"
191
"     --filter-cm          Perform common-mode noise subtraction on images\n"
Thomas White's avatar
Thomas White committed
192
193
"                           before proceeding.  Intensities will be extracted\n"
"                           from the image as it is after this processing.\n"
194
"     --filter-noise       Apply an aggressive noise filter which sets all\n"
Thomas White's avatar
Thomas White committed
195
196
197
"                           pixels in each 3x3 region to zero if any of them\n"
"                           have negative values.  Intensity measurement will\n"
"                           be performed on the image as it was before this.\n"
198
"     --unpolarized        Don't correct for the polarisation of the X-rays.\n"
199
200
"     --no-sat-corr        Don't correct values of saturated peaks using a\n"
"                           table included in the HDF5 file.\n"
201
"     --threshold=<n>      Only accept peaks above <n> ADU.  Default: 800.\n"
202
203
"     --min-gradient=<n>   Minimum gradient for Zaefferer peak search.\n"
"                           Default: 100,000.\n"
204
205
"     --min-snr=<n>        Minimum signal-to-noise ratio for peaks.\n"
"                           Default: 5.\n"
206
207
"     --min-integration-snr=<n>  Minimum signal-to-noise ratio for peaks\n"
"                                 during integration. Default: -infinity.\n"
Thomas White's avatar
Thomas White committed
208
" -e, --image=<element>    Use this image from the HDF5 file.\n"
Thomas White's avatar
Thomas White committed
209
"                           Example: /data/data0.\n"
Thomas White's avatar
Thomas White committed
210
"                           Default: The first one found.\n"
Thomas White's avatar
Thomas White committed
211
"\n"
212
213
214
215
"\nFor time-resolved stuff, you might want to use:\n\n"
"     --copy-hdf5-field <f>  Copy the value of field <f> into the stream. You\n"
"                             can use this option as many times as you need.\n"
"\n"
216
217
218
"\nOptions for greater performance or verbosity:\n\n"
"     --verbose            Be verbose about indexing.\n"
" -j <n>                   Run <n> analyses in parallel.  Default 1.\n"
219
220
221
"\n"
"\nOptions you probably won't need:\n\n"
"     --no-check-prefix    Don't attempt to correct the --prefix.\n"
222
223
224
"     --no-closer-peak     Don't integrate from the location of a nearby peak\n"
"                           instead of the position closest to the reciprocal\n"
"                           lattice point.\n"
225
226
"     --insane             Don't check that the reduced cell accounts for at\n"
"                           least 10%% of the located peaks.\n"
227
228
"     --no-bg-sub          Don't subtract local background estimates from\n"
"                           integrated intensities.\n"
229
"\n"
Thomas White's avatar
Thomas White committed
230
231
"\nYou can tune the CPU affinities for enhanced performance on NUMA machines:\n"
"\n"
Thomas White's avatar
Thomas White committed
232
233
"     --cpus=<n>           Specify number of CPUs.  This is NOT the same as\n"
"                           giving the number of analyses to run in parallel.\n"
234
235
"     --cpugroup=<n>       Batch threads in groups of this size.\n"
"     --cpuoffset=<n>      Start using CPUs at this group number.\n"
Thomas White's avatar
Thomas White committed
236
);
Thomas White's avatar
Thomas White committed
237
238
239
}


240
static void process_image(void *pp, int cookie)
241
{
242
	struct index_args *pargs = pp;
243
244
245
246
	struct hdfile *hdfile;
	struct image image;
	float *data_for_measurement;
	size_t data_size;
247
	char *filename = pargs->filename;
248
249
250
251
252
	UnitCell *cell = pargs->static_args.cell;
	int config_cmfilter = pargs->static_args.config_cmfilter;
	int config_noisefilter = pargs->static_args.config_noisefilter;
	int config_verbose = pargs->static_args.config_verbose;
	int config_polar = pargs->static_args.config_polar;
Thomas White's avatar
Thomas White committed
253
	IndexingMethod *indm = pargs->static_args.indm;
254
	struct beam_params *beam = pargs->static_args.beam;
255
256
257

	image.features = NULL;
	image.data = NULL;
Thomas White's avatar
Thomas White committed
258
	image.flags = NULL;
259
	image.indexed_cell = NULL;
260
	image.id = cookie;
261
	image.filename = filename;
262
	image.det = copy_geom(pargs->static_args.det);
263
	image.copyme = pargs->static_args.copyme;
264
265
266
267
268
269
270
	image.beam = beam;

	if ( beam == NULL ) {
		ERROR("Warning: no beam parameters file.\n");
		ERROR("I'm going to assume 1 ADU per photon, which is almost");
		ERROR(" certainly wrong.  Peak sigmas will be incorrect.\n");
	}
271

Thomas White's avatar
Thomas White committed
272
	pargs->indexable = 0;
273

274
	hdfile = hdfile_open(filename);
Thomas White's avatar
Thomas White committed
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
	if ( hdfile == NULL ) return;

	if ( pargs->static_args.element != NULL ) {

		int r;
		r = hdfile_set_image(hdfile, pargs->static_args.element);
		if ( r ) {
			ERROR("Couldn't select path '%s'\n",
			      pargs->static_args.element);
			hdfile_close(hdfile);
			return;
		}

	} else {

		int r;
		r = hdfile_set_first_image(hdfile, "/");
		if ( r ) {
			ERROR("Couldn't select first path\n");
			hdfile_close(hdfile);
			return;
		}

298
299
	}

300
	hdf5_read(hdfile, &image, pargs->static_args.config_satcorr);
301
302
303
304
305
306
307
308
309
310
311
312
313
314

	if ( (image.width != image.det->max_fs+1)
	  || (image.height != image.det->max_ss+1) )
	{
		ERROR("Image size doesn't match geometry size"
		      " - rejecting image.\n");
		ERROR("Image size: %i,%i.  Geometry size: %i,%i\n",
		      image.width, image.height,
		      image.det->max_fs+1, image.det->max_ss+1);
		hdfile_close(hdfile);
		free_detector_geometry(image.det);
		return;
	}

315
316
	if ( image.lambda < 0.0 ) {
		if ( beam != NULL ) {
Thomas White's avatar
Thomas White committed
317
			ERROR("Using nominal photon energy of %.2f eV\n",
318
                              beam->photon_energy);
Thomas White's avatar
Thomas White committed
319
320
			image.lambda = ph_en_to_lambda(
			                          eV_to_J(beam->photon_energy));
321
322
323
324
		} else {
			ERROR("No wavelength in file, so you need to give "
			      "a beam parameters file with -b.\n");
			hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
325
			free_detector_geometry(image.det);
326
327
328
			return;
		}
	}
329
	fill_in_values(image.det, hdfile);
330
331
332
333
334
335
336
337
338
339
340
341
342

	if ( config_cmfilter ) {
		filter_cm(&image);
	}

	/* Take snapshot of image after CM subtraction but before
	 * the aggressive noise filter. */
	data_size = image.width*image.height*sizeof(float);
	data_for_measurement = malloc(data_size);

	if ( config_noisefilter ) {
		filter_noise(&image, data_for_measurement);
	} else {
343
		memcpy(data_for_measurement, image.data, data_size);
344
345
	}

346
	switch ( pargs->static_args.peaks )
347
348
349
	{
	case PEAK_HDF5 :
		/* Get peaks from HDF5 */
350
351
352
		if ( get_peaks(&image, hdfile,
		               pargs->static_args.hdf5_peak_path) )
		{
353
354
355
356
			ERROR("Failed to get peaks from HDF5 file.\n");
		}
		break;
	case PEAK_ZAEF :
357
		search_peaks(&image, pargs->static_args.threshold,
358
359
		             pargs->static_args.min_gradient,
		             pargs->static_args.min_snr);
360
361
		break;
	}
Thomas White's avatar
Thomas White committed
362
363
364
365
366
367

	/* Get rid of noise-filtered version at this point
	 * - it was strictly for the purposes of peak detection. */
	free(image.data);
	image.data = data_for_measurement;

368
	/* Calculate orientation matrix (by magic) */
369
370
371
	image.div = beam->divergence;
	image.bw = beam->bandwidth;
	image.profile_radius = 0.0001e9;
372
	index_pattern(&image, cell, indm, pargs->static_args.cellr,
373
		      config_verbose, pargs->static_args.ipriv,
374
		      pargs->static_args.config_insane, pargs->static_args.tols);
375

376
	if ( image.indexed_cell != NULL ) {
Thomas White's avatar
Thomas White committed
377

378
		pargs->indexable = 1;
379

380
381
		image.reflections = find_intersections(&image,
		                                       image.indexed_cell);
382

383
		if ( image.reflections != NULL ) {
384

385
386
			integrate_reflections(&image, config_polar,
					      pargs->static_args.config_closer,
387
388
					      pargs->static_args.config_bgsub,
					      pargs->static_args.min_int_snr);
389

390
		}
391

392
393
394
395
396
	} else {

		image.reflections = NULL;

	}
Thomas White's avatar
Thomas White committed
397
398

	pthread_mutex_lock(pargs->static_args.output_mutex);
399
	write_chunk(pargs->static_args.ofh, &image, hdfile,
Thomas White's avatar
Thomas White committed
400
401
	            pargs->static_args.stream_flags);
	pthread_mutex_unlock(pargs->static_args.output_mutex);
402
403

	/* Only free cell if found */
404
	cell_free(image.indexed_cell);
405

Thomas White's avatar
Thomas White committed
406
	reflist_free(image.reflections);
407
	free(image.data);
Thomas White's avatar
Thomas White committed
408
	if ( image.flags != NULL ) free(image.flags);
409
410
	image_feature_list_free(image.features);
	hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
411
	free_detector_geometry(image.det);
412
413
414
}


415
static void *get_image(void *qp)
416
{
417
	char *line;
418
419
420
	struct index_args *pargs;
	char *rval;
	struct queue_args *qargs = qp;
421

422
	/* Initialise new task arguments */
423
424
425
	pargs = malloc(sizeof(struct index_args));
	memcpy(&pargs->static_args, &qargs->static_args,
	       sizeof(struct static_index_args));
426

427
428
429
	/* Get the next filename */
	if ( qargs->use_this_one_instead != NULL ) {

430
		line = qargs->use_this_one_instead;
431
432
433
434
		qargs->use_this_one_instead = NULL;

	} else {

435
		line = malloc(1024*sizeof(char));
436
		rval = fgets(line, 1023, qargs->fh);
Thomas White's avatar
Thomas White committed
437
438
		if ( rval == NULL ) {
			free(pargs);
Thomas White's avatar
Thomas White committed
439
			free(line);
Thomas White's avatar
Thomas White committed
440
441
			return NULL;
		}
442
443
444
		chomp(line);

	}
445

446
447
448
449
450
451
452
453
454
455
456
457
458
	if ( qargs->config_basename ) {
		char *tmp;
		tmp = safe_basename(line);
		free(line);
		line = tmp;
	}

	pargs->filename = malloc(strlen(qargs->prefix)+strlen(line)+1);

	snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line);

	free(line);

459
460
	return pargs;
}
461
462


463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
#ifdef HAVE_CLOCK_GETTIME

static time_t get_monotonic_seconds()
{
	struct timespec tp;
	clock_gettime(CLOCK_MONOTONIC, &tp);
	return tp.tv_sec;
}

#else

/* Fallback version of the above.  The time according to gettimeofday() is not
 * monotonic, so measuring intervals based on it will screw up if there's a
 * timezone change (e.g. daylight savings) while the program is running. */
static time_t get_monotonic_seconds()
{
	struct timeval tp;
	gettimeofday(&tp, NULL);
	return tp.tv_sec;
}

#endif

486
487
488
489
static void finalise_image(void *qp, void *pp)
{
	struct queue_args *qargs = qp;
	struct index_args *pargs = pp;
490
	time_t monotonic_seconds;
491

492
	qargs->n_indexable += pargs->indexable;
493
494
	qargs->n_processed++;

495
496
	monotonic_seconds = get_monotonic_seconds();
	if ( monotonic_seconds >= qargs->t_last_stats+STATS_EVERY_N_SECONDS ) {
497
498

		STATUS("%i out of %i indexed so far,"
Thomas White's avatar
Thomas White committed
499
		       " %i out of %i since the last message.\n",
500
501
		       qargs->n_indexable, qargs->n_processed,
		       qargs->n_indexable - qargs->n_indexable_last_stats,
Thomas White's avatar
Thomas White committed
502
		       qargs->n_processed - qargs->n_processed_last_stats);
503
504
505

		qargs->n_processed_last_stats = qargs->n_processed;
		qargs->n_indexable_last_stats = qargs->n_indexable;
506
		qargs->t_last_stats = monotonic_seconds;
507
508

	}
509

510
511
	free(pargs->filename);
	free(pargs);
512
513
514
}


515
516
517
static int parse_cell_reduction(const char *scellr, int *err,
                                int *reduction_needs_cell)
{
Thomas White's avatar
Thomas White committed
518
	*err = 0;
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
	if ( strcmp(scellr, "none") == 0 ) {
		*reduction_needs_cell = 0;
		return CELLR_NONE;
	} else if ( strcmp(scellr, "reduce") == 0) {
		*reduction_needs_cell = 1;
		return CELLR_REDUCE;
	} else if ( strcmp(scellr, "compare") == 0) {
		*reduction_needs_cell = 1;
		return CELLR_COMPARE;
	} else if ( strcmp(scellr, "compare_ab") == 0) {
		*reduction_needs_cell = 1;
		return CELLR_COMPARE_AB;
	} else {
		*err = 1;
		*reduction_needs_cell = 0;
		return CELLR_NONE;
	}
}


Thomas White's avatar
Thomas White committed
539
540
541
542
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
543
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
544
	FILE *fh;
545
	FILE *ofh;
546
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
547
	int n_images;
548
	int config_noindex = 0;
549
550
	int config_cmfilter = 0;
	int config_noisefilter = 0;
551
	int config_verbose = 0;
Thomas White's avatar
Thomas White committed
552
	int config_polar = 1;
553
	int config_satcorr = 1;
554
	int config_checkprefix = 1;
555
	int config_closer = 1;
556
	int config_insane = 0;
557
	int config_bgsub = 1;
558
	int config_basename = 0;
559
	float threshold = 800.0;
560
	float min_gradient = 100000.0;
Thomas White's avatar
Thomas White committed
561
	float min_snr = 5.0;
562
	double min_int_snr = -INFINITY;
Thomas White's avatar
Thomas White committed
563
564
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
565
566
567
568
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
	int indexer_needs_cell;
	int reduction_needs_cell;
Thomas White's avatar
Thomas White committed
569
	char *indm_str = NULL;
570
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
571
	char *pdb = NULL;
572
	char *prefix = NULL;
573
	char *speaks = NULL;
574
	char *scellr = NULL;
575
	char *toler = NULL;
Thomas White's avatar
Thomas White committed
576
	float tols[4] = {5.0, 5.0, 5.0, 1.5}; /* a,b,c,angles (%,%,%,deg) */
577
	int cellr;
578
	int peaks;
579
	int nthreads = 1;
580
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
581
	char *prepare_line;
582
	char prepare_filename[1024];
583
	struct queue_args qargs;
584
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
585
	char *element = NULL;
586
	double nominal_photon_energy;
Thomas White's avatar
Thomas White committed
587
	int stream_flags = STREAM_INTEGRATED;
588
589
590
591
	int cpu_num = 0;
	int cpu_groupsize = 1;
	int cpu_offset = 0;
	char *endptr;
592
	char *hdf5_peak_path = NULL;
593
594
595
596
597
598
599
	struct copy_hdf5_field *copyme;

	copyme = new_copy_hdf5_field_list();
	if ( copyme == NULL ) {
		ERROR("Couldn't allocate HDF5 field list.\n");
		return 1;
	}
Thomas White's avatar
Thomas White committed
600
601
602
603
604

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
605
		{"output",             1, NULL,               'o'},
606
		{"no-index",           0, &config_noindex,     1},
607
		{"peaks",              1, NULL,                2},
608
		{"cell-reduction",     1, NULL,                3},
609
		{"tolerance",          1, NULL,               13},
Thomas White's avatar
Thomas White committed
610
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
611
		{"geometry",           1, NULL,               'g'},
612
		{"beam",               1, NULL,               'b'},
613
614
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
615
		{"verbose",            0, &config_verbose,     1},
Thomas White's avatar
Thomas White committed
616
		{"pdb",                1, NULL,               'p'},
617
		{"prefix",             1, NULL,               'x'},
Thomas White's avatar
Thomas White committed
618
		{"unpolarized",        0, &config_polar,       0},
619
620
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
621
		{"threshold",          1, NULL,               't'},
622
		{"min-gradient",       1, NULL,                4},
623
		{"min-snr",            1, NULL,               11},
624
		{"min-integration-snr",1, NULL,               12},
625
		{"no-check-prefix",    0, &config_checkprefix, 0},
626
		{"no-closer-peak",     0, &config_closer,      0},
627
		{"insane",             0, &config_insane,      1},
Thomas White's avatar
Thomas White committed
628
		{"image",              1, NULL,               'e'},
629
		{"basename",           0, &config_basename,    1},
Thomas White's avatar
Thomas White committed
630
		{"record",             1, NULL,                5},
631
632
		{"cpus",               1, NULL,                6},
		{"cpugroup",           1, NULL,                7},
633
		{"cpuoffset",          1, NULL,                8},
Thomas White's avatar
Thomas White committed
634
		{"bg-sub",             0, &config_bgsub,       1}, /* Compat */
635
		{"no-bg-sub",          0, &config_bgsub,       0},
636
		{"hdf5-peaks",         1, NULL,                9},
637
		{"copy-hdf5-field",    1, NULL,               10},
Thomas White's avatar
Thomas White committed
638
639
640
641
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
642
	while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:b:e:",
Thomas White's avatar
Thomas White committed
643
	                        longopts, NULL)) != -1) {
Thomas White's avatar
Thomas White committed
644
645

		switch (c) {
Thomas White's avatar
Thomas White committed
646
		case 'h' :
Thomas White's avatar
Thomas White committed
647
648
649
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
650
		case 'i' :
Thomas White's avatar
Thomas White committed
651
652
653
			filename = strdup(optarg);
			break;

654
655
656
657
		case 'o' :
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
658
		case 'z' :
Thomas White's avatar
Thomas White committed
659
660
661
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
662
		case 'p' :
Thomas White's avatar
Thomas White committed
663
664
665
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
666
		case 'x' :
667
668
669
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
670
		case 'j' :
671
672
673
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
674
		case 'g' :
Thomas White's avatar
Thomas White committed
675
676
677
			geometry = strdup(optarg);
			break;

678
679
680
681
		case 't' :
			threshold = strtof(optarg, NULL);
			break;

682
683
684
685
686
687
688
689
690
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

691
692
693
694
		case 2 :
			speaks = strdup(optarg);
			break;

695
696
697
698
		case 3 :
			scellr = strdup(optarg);
			break;

699
700
701
702
		case 13 :
			toler = strdup(optarg);
			break;

703
704
705
706
		case 4 :
			min_gradient = strtof(optarg, NULL);
			break;

707
708
709
710
		case 11 :
			min_snr = strtof(optarg, NULL);
			break;

711
712
713
714
		case 12 :
			min_int_snr = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
715
716
717
718
		case 'e' :
			element = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
719
720
721
722
723
		case 5 :
			stream_flags = parse_stream_flags(optarg);
			if ( stream_flags < 0 ) return 1;
			break;

724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
		case 6 :
			cpu_num = strtol(optarg, &endptr, 10);
			if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
				ERROR("Invalid number of CPUs ('%s')\n",
				      optarg);
				return 1;
			}
			break;

		case 7 :
			cpu_groupsize = strtol(optarg, &endptr, 10);
			if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
				ERROR("Invalid CPU group size ('%s')\n",
				      optarg);
				return 1;
			}
			if ( cpu_groupsize < 1 ) {
				ERROR("CPU group size cannot be"
				      " less than 1.\n");
				return 1;
			}
			break;

		case 8 :
			cpu_offset = strtol(optarg, &endptr, 10);
			if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
				ERROR("Invalid CPU offset ('%s')\n",
				      optarg);
				return 1;
			}
			if ( cpu_offset < 0 ) {
				ERROR("CPU offset must be positive.\n");
				return 1;
			}
			break;

760
761
762
763
		case 9 :
			hdf5_peak_path = strdup(optarg);
			break;

764
765
766
767
		case 10 :
			add_copy_hdf5_field(copyme, optarg);
			break;

Thomas White's avatar
Thomas White committed
768
		case 0 :
Thomas White's avatar
Thomas White committed
769
770
			break;

Thomas White's avatar
Thomas White committed
771
		default :
Thomas White's avatar
Thomas White committed
772
773
774
775
776
			return 1;
		}

	}

777
778
779
780
781
782
783
	if ( (cpu_num > 0) && (cpu_num % cpu_groupsize != 0) ) {
		ERROR("Number of CPUs must be divisible by"
		      " the CPU group size.\n");
		return 1;
	}


Thomas White's avatar
Thomas White committed
784
785
786
787
788
789
790
791
792
	if ( filename == NULL ) {
		filename = strdup("-");
	}
	if ( strcmp(filename, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(filename, "r");
	}
	if ( fh == NULL ) {
Thomas White's avatar
Thomas White committed
793
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
794
795
		return 1;
	}
Thomas White's avatar
Thomas White committed
796
	free(filename);
Thomas White's avatar
Thomas White committed
797

798
799
800
801
802
803
804
805
806
807
808
809
810
811
	if ( outfile == NULL ) {
		outfile = strdup("-");
	}
	if ( strcmp(outfile, "-") == 0 ) {
		ofh = stdout;
	} else {
		ofh = fopen(outfile, "w");
	}
	if ( ofh == NULL ) {
		ERROR("Failed to open output file '%s'\n", outfile);
		return 1;
	}
	free(outfile);

812
813
814
815
	if ( hdf5_peak_path == NULL ) {
		hdf5_peak_path = strdup("/processing/hitfinder/peakinfo");
	}

816
817
818
819
820
821
822
823
824
825
826
827
828
	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;
	}
Thomas White's avatar
Thomas White committed
829
	free(speaks);
830

Thomas White's avatar
Thomas White committed
831
832
833
834
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

835
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
836
		prefix = strdup("");
837
	} else {
838
839
840
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
841
842
	}

843
	if ( nthreads == 0 ) {
844
845
846
847
		ERROR("Invalid number of threads.\n");
		return 1;
	}

848
849
	if ( (indm_str == NULL) ||
	     ((indm_str != NULL) && (strcmp(indm_str, "none") == 0)) ) {
850
		STATUS("Not indexing anything.\n");
851
		indexer_needs_cell = 0;
852
853
854
		reduction_needs_cell = 0;
		indm = NULL;
		cellr = CELLR_NONE;
Thomas White's avatar
Thomas White committed
855
	} else {
856
857
858
859
860
861
862
863
		if ( indm_str == NULL ) {
			STATUS("You didn't specify an indexing method, so I "
			       " won't try to index anything.\n"
			       "If that isn't what you wanted, re-run with"
			       " --indexing=<method>.\n");
			indm = NULL;
			indexer_needs_cell = 0;
		} else {
Thomas White's avatar
Thomas White committed
864
865
			indm = build_indexer_list(indm_str,
			                          &indexer_needs_cell);
866
867
868
869
870
			if ( indm == NULL ) {
				ERROR("Invalid indexer list '%s'\n", indm_str);
				return 1;
			}
			free(indm_str);
Thomas White's avatar
Thomas White committed
871
		}
Thomas White's avatar
Thomas White committed
872

873
874
875
876
877
878
879
		reduction_needs_cell = 0;
		if ( scellr == NULL ) {
			STATUS("You didn't specify a cell reduction method, so"
			       " I'm going to use 'reduce'.\n");
			cellr = CELLR_REDUCE;
			reduction_needs_cell = 1;
		} else {
880
881
882
883
884
885
886
887
888
			int err;
			cellr = parse_cell_reduction(scellr, &err,
			                             &reduction_needs_cell);
			if ( err ) {
				ERROR("Unrecognised cell reduction '%s'\n",
			              scellr);
				return 1;
			}
			free(scellr);
889
		}
890
891
	}

892
893
894
	/* No indexing -> no reduction */
	if ( indm == NULL ) reduction_needs_cell = 0;

895
896
	if ( toler != NULL ) {
		int ttt;
Thomas White's avatar
Thomas White committed
897
898
		ttt = sscanf(toler, "%f,%f,%f,%f",
		             &tols[0], &tols[1], &tols[2], &tols[3] );
899
900
901
902
903
904
		if ( ttt != 4 ) {
			ERROR("Invalid parameters for '--tolerance'\n");
			return 1;
		}
	}

Thomas White's avatar
Thomas White committed
905
906
907
908
909
910
911
912
913
914
915
916
	if ( geometry == NULL ) {
		ERROR("You need to specify a geometry file with --geometry\n");
		return 1;
	}

	det = get_detector_geometry(geometry);
	if ( det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geometry);
		return 1;
	}
	free(geometry);

Thomas White's avatar
Thomas White committed
917
	if ( reduction_needs_cell || indexer_needs_cell ) {
918
919
		cell = load_cell_from_pdb(pdb);
		if ( cell == NULL ) {
920
921
922
			ERROR("Couldn't read unit cell (from %s)\n", pdb);
			return 1;
		}
923
	} else {
Thomas White's avatar
Thomas White committed
924
925
		STATUS("No cell needed for these choices of indexing"
		       " and reduction.\n");
926
		cell = NULL;
927
	}
928
	free(pdb);
929

Thomas White's avatar
Thomas White committed
930
	write_stream_header(ofh, argc, argv);
931

932
933
934
935
936
937
938
939
	if ( beam != NULL ) {
		nominal_photon_energy = beam->photon_energy;
	} else {
		STATUS("No beam parameters file was given, so I'm taking the"
		       " nominal photon energy to be 2 keV.\n");
		nominal_photon_energy = 2000.0;
	}

940
	/* Get first filename and use it to set up the indexing */
941
	prepare_line = malloc(1024*sizeof(char));
942
943
944
945
	rval = fgets(prepare_line, 1023, fh);
	if ( rval == NULL ) {
		ERROR("Failed to get filename to prepare indexing.\n");
		return 1;
946
	}
947
	chomp(prepare_line);
948
949
950
951
952
953
	if ( config_basename ) {
		char *tmp;
		tmp = safe_basename(prepare_line);
		free(prepare_line);
		prepare_line = tmp;
	}
954
955
956
957
	snprintf(prepare_filename, 1023, "%s%s", prefix, prepare_line);
	qargs.use_this_one_instead = prepare_line;

	/* Prepare the indexer */
Thomas White's avatar
Thomas White committed
958
959
960
961
962
963
964
965
966
	if ( indm != NULL ) {
		ipriv = prepare_indexing(indm, cell, prepare_filename, det,
		                         nominal_photon_energy);
		if ( ipriv == NULL ) {
			ERROR("Failed to prepare indexing.\n");
			return 1;
		}
	} else {
		ipriv = NULL;
967
968
	}

Thomas White's avatar
Thomas White committed
969
	gsl_set_error_handler_off();
Thomas White's avatar
Thomas White committed
970

971
972
973
974
975
976
977
	qargs.static_args.cell = cell;
	qargs.static_args.config_cmfilter = config_cmfilter;
	qargs.static_args.config_noisefilter = config_noisefilter;
	qargs.static_args.config_verbose = config_verbose;
	qargs.static_args.config_polar = config_polar;
	qargs.static_args.config_satcorr = config_satcorr;
	qargs.static_args.config_closer = config_closer;
978
	qargs.static_args.config_insane = config_insane;
979
	qargs.static_args.config_bgsub = config_bgsub;
980
	qargs.static_args.cellr = cellr;
981
982
983
984
	qargs.static_args.tols[0] = tols[0];
	qargs.static_args.tols[1] = tols[1];
	qargs.static_args.tols[2] = tols[2];
	qargs.static_args.tols[3] = tols[3];
985
	qargs.static_args.threshold = threshold;
986
	qargs.static_args.min_gradient = min_gradient;
987
	qargs.static_args.min_snr = min_snr;
Thomas White's avatar
Thomas White committed
988
	qargs.static_args.min_int_snr = min_int_snr;
989
990
991
992
993
994
	qargs.static_args.det = det;
	qargs.static_args.indm = indm;
	qargs.static_args.ipriv = ipriv;
	qargs.static_args.peaks = peaks;
	qargs.static_args.output_mutex = &output_mutex;
	qargs.static_args.ofh = ofh;
995
	qargs.static_args.beam = beam;
Thomas White's avatar
Thomas White committed
996
	qargs.static_args.element = element;
997
	qargs.static_args.stream_flags = stream_flags;
998
	qargs.static_args.hdf5_peak_path = hdf5_peak_path;
999
	qargs.static_args.copyme = copyme;
1000

For faster browsing, not all history is shown. View entire blame