indexamajig.c 34.7 KB
Newer Older
Thomas White's avatar
Thomas White committed
1
/*
Chuck's avatar
Chuck committed
2
 * indexamajigFork.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
 *
Chuck's avatar
Chuck committed
6
7
 * (c) 2006-2012 Thomas White <taw@physics.org>
 * (c) 2012-     Chunhong Yoon <chun.hong.yoon@desy.de>
8
 *
Chuck's avatar
Chuck committed
9
 * Part of CrystFEL - crystallography with a FEL
Thomas White's avatar
Thomas White committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 *
 */


#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
24
#include <hdf5.h>
Thomas White's avatar
Thomas White committed
25
#include <gsl/gsl_errno.h>
26
#include <pthread.h>
Chuck's avatar
Chuck committed
27
28
#include <sys/wait.h>
#include <fcntl.h>
29
30

#ifdef HAVE_CLOCK_GETTIME
31
#include <time.h>
32
33
34
#else
#include <sys/time.h>
#endif
Thomas White's avatar
Thomas White committed
35

Chuck's avatar
Chuck committed
36
37
38
39
40
41
42
43
44
45
46
#include <crystfel/utils.h>
#include <crystfel/hdf5-file.h>
#include <crystfel/index.h>
#include <crystfel/peaks.h>
#include <crystfel/detector.h>
#include <crystfel/filters.h>
#include <crystfel/thread-pool.h>
#include <crystfel/beam-parameters.h>
#include <crystfel/geometry.h>
#include <crystfel/stream.h>
#include <crystfel/reflist-utils.h>
47

48

Thomas White's avatar
Thomas White committed
49
/* Write statistics at APPROXIMATELY this interval */
Chuck's avatar
Chuck committed
50
51
52
#define STATS_EVERY_N_SECONDS (2)

#define LINE_LENGTH 1024
53

Chuck's avatar
Chuck committed
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#define BUFSIZE 64

#ifdef HAVE_CLOCK_GETTIME
static double get_time()
{
	struct timespec tp;
	clock_gettime(CLOCK_MONOTONIC, &tp);
	double sec = (double) tp.tv_sec+ (double) tp.tv_nsec/1000000000;
	return sec; //nano resolution
}
#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 double get_time()
{
	struct timeval tp;
	gettimeofday(&tp, NULL);
	double sec = (double) tp.tv_sec+ (double) tp.tv_usec/1000000;
	return sec; //micro resolution
}
#endif
76

77
78
79
80
81
82
enum {
	PEAK_ZAEF,
	PEAK_HDF5,
};


83
84
/* Information about the indexing process which is common to all patterns */
struct static_index_args
85
86
87
88
89
{
	UnitCell *cell;
	int config_cmfilter;
	int config_noisefilter;
	int config_verbose;
Thomas White's avatar
Thomas White committed
90
	int stream_flags;         /* What goes into the output? */
Chuck's avatar
Chuck committed
91
	int config_polar;
92
	int config_satcorr;
93
	int config_closer;
94
	int config_insane;
95
	int config_bgsub;
96
	float threshold;
97
	float min_gradient;
98
	float min_snr;
99
	double min_int_snr;
Thomas White's avatar
Thomas White committed
100
	struct detector *det;
Thomas White's avatar
Thomas White committed
101
102
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
Thomas White's avatar
Thomas White committed
103
	int peaks;                /* Peak detection method */
104
	int cellr;
105
	float tols[4];
106
	struct beam_params *beam;
Thomas White's avatar
Thomas White committed
107
	const char *element;
108
	const char *hdf5_peak_path;
109
110
111
	double ir_inn;
	double ir_mid;
	double ir_out;
Chuck's avatar
Chuck committed
112
	
113
114
115
	/* Output stream */
	pthread_mutex_t *output_mutex;  /* Protects the output stream */
	FILE *ofh;
Chuck's avatar
Chuck committed
116
117

	char *outfile;	
118
	const struct copy_hdf5_field *copyme;
Chuck's avatar
Chuck committed
119
	int nProcesses;
120
121
122
};


123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
/* 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;
140
	int config_basename;
141
142
	struct static_index_args static_args;

143
	char *use_this_one_instead;
144
145
146
147
148

	int n_indexable;
	int n_processed;
	int n_indexable_last_stats;
	int n_processed_last_stats;
Chuck's avatar
Chuck committed
149
150
151
152
153
154
155
156
157
158
159

	int n_indexableTotal;
	int n_processedTotal;
	int n_indexable_last_statsTotal;
	int n_processed_last_statsTotal;

	int nPerProcess;
	int updateReader;

	double t_last_stats;

160
161
};

Chuck's avatar
Chuck committed
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
// Count number of patterns in .lst
int count_patterns(FILE *fh){
    char *rval;
    int n_total = 0;
    do{
        char line[LINE_LENGTH];
        rval = fgets(line,LINE_LENGTH-1,fh);
        if (rval != NULL){
            n_total++;
        }
    }while(rval!=NULL);
    if (ferror(fh)) {
        printf("Read error\n");
        return -1;
    }
    return n_total;
}

// Assign a batch number to a process
int getBatchNum(int pid[], int length){
    int i, id = 0;
    for (i=0; i<length; i++) {
        id += ((pid[i]>0)*pow(2,i)); // child = 1, parent = 0
    }
    return id;
}
188

Thomas White's avatar
Thomas White committed
189
190
191
192
193
194
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
"Process and index FEL diffraction images.\n"
"\n"
195
" -h, --help               Display this help message.\n"
Thomas White's avatar
Thomas White committed
196
"\n"
197
" -i, --input=<filename>   Specify file containing list of images to process.\n"
Thomas White's avatar
Thomas White committed
198
"                           '-' means stdin, which is the default.\n"
Thomas White's avatar
Thomas White committed
199
200
" -o, --output=<filename>  Write output stream to this file. '-' for stdout.\n"
"                           Default: indexamajig.stream\n"
Thomas White's avatar
Thomas White committed
201
"\n"
Thomas White's avatar
Thomas White committed
202
203
204
205
206
"     --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
207
"                            reax     : DPS algorithm with known unit cell\n"
208
" -g. --geometry=<file>    Get detector geometry from file.\n"
209
210
211
" -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"
212
" -p, --pdb=<file>         PDB file from which to get the unit cell to match.\n"
Thomas White's avatar
Thomas White committed
213
"                           Default: 'molecule.pdb'.\n"
214
"     --basename           Remove the directory parts of the filenames.\n"
215
" -x, --prefix=<p>         Prefix filenames from input file with <p>.\n"
216
217
218
"     --peaks=<method>     Use 'method' for finding peaks.  Choose from:\n"
"                           zaef  : Use Zaefferer (2000) gradient detection.\n"
"                                    This is the default method.\n"
219
220
221
"                           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
222
223
"\n\n"
"You can control what information is included in the output stream using\n"
Thomas White's avatar
Thomas White committed
224
"' --record=<flag1>,<flag2>,<flag3>' and so on.  Possible flags are:\n\n"
225
226
" integrated        Include a list of reflection intensities, produced by\n"
"                    integrating around predicted peak locations.\n"
Thomas White's avatar
Thomas White committed
227
"\n"
228
229
" peaks             Include peak locations and intensities from the peak\n"
"                    search.\n"
Thomas White's avatar
Thomas White committed
230
"\n"
231
" peaksifindexed    As 'peaks', but only if the pattern could be indexed.\n"
Thomas White's avatar
Thomas White committed
232
"\n"
233
234
" peaksifnotindexed As 'peaks', but only if the pattern could NOT be indexed.\n"
"\n\n"
Thomas White's avatar
Thomas White committed
235
"The default is '--record=integrated'.\n"
Thomas White's avatar
Thomas White committed
236
237
"\n\n"
"For more control over the process, you might need:\n\n"
Chuck's avatar
Chuck committed
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
"     --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"
"                           compare_ab : compare 'a' and 'b' lengths only.\n"
"     --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"
"     --filter-cm          Perform common-mode noise subtraction on images\n"
"                           before proceeding.  Intensities will be extracted\n"
"                           from the image as it is after this processing.\n"
"     --filter-noise       Apply an aggressive noise filter which sets all\n"
"                           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"
"     --unpolarized        Don't correct for the polarisation of the X-rays.\n"
"     --no-sat-corr        Don't correct values of saturated peaks using a\n"
"                           table included in the HDF5 file.\n"
"     --threshold=<n>      Only accept peaks above <n> ADU.  Default: 800.\n"
"     --min-gradient=<n>   Minimum gradient for Zaefferer peak search.\n"
"                           Default: 100,000.\n"
"     --min-snr=<n>        Minimum signal-to-noise ratio for peaks.\n"
"                           Default: 5.\n"
"     --min-integration-snr=<n>  Minimum signal-to-noise ratio for peaks\n"
"                                 during integration. Default: -infinity.\n"
" -e, --image=<element>    Use this image from the HDF5 file.\n"
"                           Example: /data/data0.\n"
"                           Default: The first one found.\n"
Thomas White's avatar
Thomas White committed
267
"\n"
268
269
270
271
"\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"
272
273
274
"\nOptions for greater performance or verbosity:\n\n"
"     --verbose            Be verbose about indexing.\n"
" -j <n>                   Run <n> analyses in parallel.  Default 1.\n"
275
276
277
"\n"
"\nOptions you probably won't need:\n\n"
"     --no-check-prefix    Don't attempt to correct the --prefix.\n"
278
279
280
"     --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"
281
282
"     --insane             Don't check that the reduced cell accounts for at\n"
"                           least 10%% of the located peaks.\n"
283
284
"     --no-bg-sub          Don't subtract local background estimates from\n"
"                           integrated intensities.\n"
285
"\n"
Thomas White's avatar
Thomas White committed
286
287
"\nYou can tune the CPU affinities for enhanced performance on NUMA machines:\n"
"\n"
Thomas White's avatar
Thomas White committed
288
289
"     --cpus=<n>           Specify number of CPUs.  This is NOT the same as\n"
"                           giving the number of analyses to run in parallel.\n"
290
291
"     --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
292
);
Thomas White's avatar
Thomas White committed
293
294
}

Chuck's avatar
Chuck committed
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
int readUpdate(void *qp, int fd_pipe[][2], int finish){
	struct queue_args *qargs = qp;
	int i,j;
	int rFlag;
	int pipesOpen = 0;
	char bufferR[BUFSIZE];			
	memset(bufferR, 0, BUFSIZE);
	int numFields = 2;

	for (i=0; i<qargs->static_args.nProcesses; i++){
		if (finish && i == qargs->updateReader) {
			// do nothing
		}else{
			rFlag = read(fd_pipe[i][0], bufferR, BUFSIZE-1); // wait till something to read
			if (rFlag != 0){
				char * pch; 
				char delims[] = "#";
  				pch = strtok (bufferR, delims); 
  				for (j=0;j<numFields;j++) 
  				{ 
					switch (j){
						case 0:
							qargs->n_indexableTotal += atoi(pch);	
							break;
						case 1:
							qargs->n_processedTotal += atoi(pch);	
						break;
					}
    				pch = strtok (NULL, delims); 
  				} 
				memset(bufferR, 0, BUFSIZE);
				pipesOpen++;
			}
		}
	}
	STATUS("%i out of %i indexed so far,"
		       " %i out of %i since the last message.\n\n",
			qargs->n_indexableTotal, qargs->n_processedTotal,
			qargs->n_indexableTotal-qargs->n_indexable_last_statsTotal,
			qargs->n_processedTotal-qargs->n_processed_last_statsTotal);
Thomas White's avatar
Thomas White committed
335

Chuck's avatar
Chuck committed
336
337
338
339
340
341
342
343
344
345
	qargs->n_indexable_last_statsTotal = qargs->n_indexableTotal;
	qargs->n_processed_last_statsTotal = qargs->n_processedTotal;

	return pipesOpen;	
}

// Manipulate image
static void process_image(char ***array, int batch, void *qp, int fd_pipe[][2]) {
	
	struct queue_args *qargs = qp;
346
347
	float *data_for_measurement;
	size_t data_size;
Chuck's avatar
Chuck committed
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
	UnitCell *cell = qargs->static_args.cell;
	int config_cmfilter = qargs->static_args.config_cmfilter;
	int config_noisefilter = qargs->static_args.config_noisefilter;
	int config_verbose = qargs->static_args.config_verbose;
	//int config_polar = qargs->static_args.config_polar;
	IndexingMethod *indm = qargs->static_args.indm;
	struct beam_params *beam = qargs->static_args.beam;

	int i;	
	int r, check;
	struct hdfile *hdfile;
	char *outfile = qargs->static_args.outfile;
	int nPerBatch = qargs->nPerProcess;

	for (i=0; i<nPerBatch; i++) {
		struct image image;
		image.features = NULL;
		image.data = NULL;
		image.flags = NULL;
		image.indexed_cell = NULL;
		image.det = copy_geom(qargs->static_args.det);
		image.copyme = qargs->static_args.copyme;
		image.beam = beam;
		image.id = batch;	// MUST SET ID FOR MOSFLM TO WORK PROPERLY

		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");
		}
378

Chuck's avatar
Chuck committed
379
		//pargs->indexable = 0;
Thomas White's avatar
Thomas White committed
380

Chuck's avatar
Chuck committed
381
382
383
384
385
386
387
388
389
390
391
		char *filename = NULL;
		char *line = array[batch][i];
		chomp(line);
//printf("%d ***************%s\n",batch,line);
		filename = malloc(strlen(qargs->prefix)+strlen(line)+1);
		snprintf(filename,LINE_LENGTH-1,"%s%s",qargs->prefix,line); //maximum print length
		free(line);
		image.filename = filename;
//printf("%d ***************%s\n",batch,filename);
		hdfile = hdfile_open(filename);
		if (hdfile == NULL) return;
Thomas White's avatar
Thomas White committed
392

Chuck's avatar
Chuck committed
393
394
395
		r = hdfile_set_first_image(hdfile, "/"); // Need this to read hdf5 files
		if (r) {
			ERROR("Couldn't select first path\n");
Thomas White's avatar
Thomas White committed
396
397
398
399
			hdfile_close(hdfile);
			return;
		}

Chuck's avatar
Chuck committed
400
401
		check = hdf5_read(hdfile,&image,1);
		if (check == 1){
Thomas White's avatar
Thomas White committed
402
403
404
405
			hdfile_close(hdfile);
			return;
		}

Chuck's avatar
Chuck committed
406
407
408
409
410
411
412
413
		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);
414
			hdfile_close(hdfile);
Thomas White's avatar
Thomas White committed
415
			free_detector_geometry(image.det);
416
417
			return;
		}
418

Chuck's avatar
Chuck committed
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
		if ( image.lambda < 0.0 ) {
			if ( beam != NULL ) {
				ERROR("Using nominal photon energy of %.2f eV\n",
    	                          beam->photon_energy);
				image.lambda = ph_en_to_lambda(
			                          eV_to_J(beam->photon_energy));
			} else {
				ERROR("No wavelength in file, so you need to give "
			      "a beam parameters file with -b.\n");
				hdfile_close(hdfile);
				free_detector_geometry(image.det);
				return;
			}
		}
		fill_in_values(image.det, hdfile);
434

Chuck's avatar
Chuck committed
435
436
437
		if ( config_cmfilter ) {
			filter_cm(&image);
		}
438

Chuck's avatar
Chuck committed
439
440
441
442
		/* 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);
443

Chuck's avatar
Chuck committed
444
445
446
447
		if ( config_noisefilter ) {
			filter_noise(&image, data_for_measurement);
		} else {
			memcpy(data_for_measurement, image.data, data_size);
448
		}
Thomas White's avatar
Thomas White committed
449

Chuck's avatar
Chuck committed
450
451
452
453
454
455
456
457
458
459
		switch ( qargs->static_args.peaks )
		{
		case PEAK_HDF5 :
			/* Get peaks from HDF5 */
			if ( get_peaks(&image, hdfile,
		               qargs->static_args.hdf5_peak_path) )
			{
				ERROR("Failed to get peaks from HDF5 file.\n");
			}
			break;
Thomas White's avatar
Thomas White committed
460
		case PEAK_ZAEF :
Chuck's avatar
Chuck committed
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
486
487
			search_peaks(&image, qargs->static_args.threshold,
		             qargs->static_args.min_gradient,
		             qargs->static_args.min_snr,
		             qargs->static_args.ir_inn,
		             qargs->static_args.ir_mid,
		             qargs->static_args.ir_out);
			break;
		}
		
		/* 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;

		/* Calculate orientation matrix (by magic) */
		image.div = beam->divergence;
		image.bw = beam->bandwidth;
		image.profile_radius = 0.0001e9;

		///// RUN INDEXING HERE /////
		index_pattern(&image, cell, indm, qargs->static_args.cellr,
		     config_verbose, qargs->static_args.ipriv,
		     qargs->static_args.config_insane, qargs->static_args.tols);

		if ( image.indexed_cell != NULL ) {
			//pargs->indexable = 1;
			image.reflections = find_intersections(&image,
488
		                                       image.indexed_cell);
Chuck's avatar
Chuck committed
489
490
491
492
493
494
495
496
497
498
499
			if ( image.reflections != NULL ) {
				integrate_reflections(&image,
					      qargs->static_args.config_closer,
					      qargs->static_args.config_bgsub,
					      qargs->static_args.min_int_snr,
					      qargs->static_args.ir_inn,
					      qargs->static_args.ir_mid,
					      qargs->static_args.ir_out);
			}
		} else {
			image.reflections = NULL;
500
		}
501

Chuck's avatar
Chuck committed
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
		// Write Lock
		struct flock fl = {F_WRLCK, SEEK_SET,   0,      0,     0 };
    	int fd;
	    fl.l_pid = getpid();
	
		char *outfilename = NULL;
		chomp(outfile); // may not need this
//printf("%d ***************%s\n",batch,outfile);		
		outfilename = malloc(strlen(outfile)+1);
		snprintf(outfilename,LINE_LENGTH-1,"%s",outfile); //maximum print length
//printf("%d ***************%s\n",batch,outfilename);
   		if ((fd = open(outfilename, O_WRONLY)) == -1) {
        	perror("Error on opening\n");
        	exit(1);
    	}

	    if (fcntl(fd, F_SETLKW, &fl) == -1) {
	        perror("Error on setting lock wait\n");
	        exit(1);
	    }

		// LOCKED! Write chunk
		FILE *fh;
//printf("%d ***************%s\n",batch,outfilename);
		fh = fopen(outfilename,"a");
		if (fh == NULL) { 
			perror("Error inside lock\n"); 
		}
		write_chunk(fh, &image, hdfile, qargs->static_args.stream_flags);		
		fclose(fh);

		// Unlock stream for other processes
	    fl.l_type = F_UNLCK;  // set to unlock same region
  	  	if (fcntl(fd, F_SETLK, &fl) == -1) {
   	    	perror("fcntl");
   	     	exit(1);
    	}
    	close(fd);

		///// WRITE UPDATE /////
		double seconds;
		qargs->n_indexable += ( image.indexed_cell != NULL );
		qargs->n_processed++;
		seconds = get_time();
		if ( seconds >= qargs->t_last_stats+STATS_EVERY_N_SECONDS 
			|| qargs->n_processed == qargs->nPerProcess) { // Must write if finished

			// WRITE PIPE HERE
			char bufferW[BUFSIZE];
			memset(bufferW, 0, BUFSIZE);
			sprintf(bufferW,"%d#%d",
				qargs->n_indexable - qargs->n_indexable_last_stats, 
				qargs->n_processed - qargs->n_processed_last_stats);
			write(fd_pipe[batch][1], bufferW, BUFSIZE-1);

			// Update stats
			qargs->n_processed_last_stats = qargs->n_processed;
			qargs->n_indexable_last_stats = qargs->n_indexable;
			qargs->t_last_stats = seconds;

			///// READ UPDATE /////
			if (batch == qargs->updateReader) {
				readUpdate(qargs,fd_pipe, 0);
			}
Thomas White's avatar
Thomas White committed
566
		}
567

Chuck's avatar
Chuck committed
568
569
570
571
572
573
574
575
576
		///// FREE /////
		/* Only free cell if found */
		cell_free(image.indexed_cell);
		reflist_free(image.reflections);
		free(image.data);
		if ( image.flags != NULL ) free(image.flags);
		image_feature_list_free(image.features);
		hdfile_close(hdfile);
		free_detector_geometry(image.det);
577
	}
578
579
}

580
581
582
static int parse_cell_reduction(const char *scellr, int *err,
                                int *reduction_needs_cell)
{
Thomas White's avatar
Thomas White committed
583
	*err = 0;
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
	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
604
605
606
607
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
608
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
609
	FILE *fh;
610
	FILE *ofh;
611
	char *rval = NULL;
Thomas White's avatar
Thomas White committed
612
	int n_images;
613
	int config_noindex = 0;
614
615
	int config_cmfilter = 0;
	int config_noisefilter = 0;
616
	int config_verbose = 0;
Chuck's avatar
Chuck committed
617
	int config_polar = 1;
618
	int config_satcorr = 1;
619
	int config_checkprefix = 1;
620
	int config_closer = 1;
621
	int config_insane = 0;
622
	int config_bgsub = 1;
623
	int config_basename = 0;
624
	float threshold = 800.0;
625
	float min_gradient = 100000.0;
Thomas White's avatar
Thomas White committed
626
	float min_snr = 5.0;
627
	double min_int_snr = -INFINITY;
Thomas White's avatar
Thomas White committed
628
629
	struct detector *det;
	char *geometry = NULL;
Thomas White's avatar
Thomas White committed
630
631
632
633
	IndexingMethod *indm;
	IndexingPrivate **ipriv;
	int indexer_needs_cell;
	int reduction_needs_cell;
Thomas White's avatar
Thomas White committed
634
	char *indm_str = NULL;
635
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
636
	char *pdb = NULL;
637
	char *prefix = NULL;
638
	char *speaks = NULL;
639
	char *scellr = NULL;
640
	char *toler = NULL;
Thomas White's avatar
Thomas White committed
641
	float tols[4] = {5.0, 5.0, 5.0, 1.5}; /* a,b,c,angles (%,%,%,deg) */
642
	int cellr;
643
	int peaks;
Chuck's avatar
Chuck committed
644
	int nProcesses = 1;
645
	pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
646
	char *prepare_line;
647
	char prepare_filename[1024];
648
	struct queue_args qargs;
649
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
650
	char *element = NULL;
651
	double nominal_photon_energy;
Thomas White's avatar
Thomas White committed
652
	int stream_flags = STREAM_INTEGRATED;
653
654
655
656
	int cpu_num = 0;
	int cpu_groupsize = 1;
	int cpu_offset = 0;
	char *endptr;
657
	char *hdf5_peak_path = NULL;
658
	struct copy_hdf5_field *copyme;
659
660
661
662
	char *intrad = NULL;
	float ir_inn = 4.0;
	float ir_mid = 5.0;
	float ir_out = 7.0;
Chuck's avatar
Chuck committed
663
	
664
665
666
667
668
	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
669
670
671
672
673

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
674
		{"output",             1, NULL,               'o'},
675
		{"no-index",           0, &config_noindex,     1},
Thomas White's avatar
Thomas White committed
676
		{"indexing",           1, NULL,               'z'},
Thomas White's avatar
Thomas White committed
677
		{"geometry",           1, NULL,               'g'},
678
		{"beam",               1, NULL,               'b'},
679
680
		{"filter-cm",          0, &config_cmfilter,    1},
		{"filter-noise",       0, &config_noisefilter, 1},
681
		{"verbose",            0, &config_verbose,     1},
Thomas White's avatar
Thomas White committed
682
		{"pdb",                1, NULL,               'p'},
683
		{"prefix",             1, NULL,               'x'},
684
685
		{"no-sat-corr",        0, &config_satcorr,     0},
		{"sat-corr",           0, &config_satcorr,     1}, /* Compat */
686
		{"threshold",          1, NULL,               't'},
687
		{"no-check-prefix",    0, &config_checkprefix, 0},
688
		{"no-closer-peak",     0, &config_closer,      0},
689
		{"insane",             0, &config_insane,      1},
Thomas White's avatar
Thomas White committed
690
		{"image",              1, NULL,               'e'},
691
		{"basename",           0, &config_basename,    1},
Thomas White's avatar
Thomas White committed
692
693
694
695
696
697
		{"bg-sub",             0, &config_bgsub,       1}, /* Compat */
		{"no-bg-sub",          0, &config_bgsub,       0},

		{"peaks",              1, NULL,                2},
		{"cell-reduction",     1, NULL,                3},
		{"min-gradient",       1, NULL,                4},
Thomas White's avatar
Thomas White committed
698
		{"record",             1, NULL,                5},
699
700
		{"cpus",               1, NULL,                6},
		{"cpugroup",           1, NULL,                7},
701
		{"cpuoffset",          1, NULL,                8},
702
		{"hdf5-peaks",         1, NULL,                9},
703
		{"copy-hdf5-field",    1, NULL,               10},
Thomas White's avatar
Thomas White committed
704
705
706
		{"min-snr",            1, NULL,               11},
		{"min-integration-snr",1, NULL,               12},
		{"tolerance",          1, NULL,               13},
707
		{"int-radius",         1, NULL,               14},
Thomas White's avatar
Thomas White committed
708
709
710
711
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
712
	while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:b:e:",
Thomas White's avatar
Thomas White committed
713
714
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
715
		switch (c) {
Thomas White's avatar
Thomas White committed
716
717

			case 'h' :
Thomas White's avatar
Thomas White committed
718
719
720
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
721
			case 'i' :
Thomas White's avatar
Thomas White committed
722
723
724
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
725
			case 'o' :
726
727
728
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
729
			case 'z' :
Thomas White's avatar
Thomas White committed
730
731
732
			indm_str = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
733
			case 'p' :
Thomas White's avatar
Thomas White committed
734
735
736
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
737
			case 'x' :
738
739
740
			prefix = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
741
			case 'j' :
Chuck's avatar
Chuck committed
742
			nProcesses = atoi(optarg);
743
744
			break;

Thomas White's avatar
Thomas White committed
745
			case 'g' :
Thomas White's avatar
Thomas White committed
746
747
748
			geometry = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
749
			case 't' :
750
751
752
			threshold = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
753
			case 'b' :
754
755
756
757
758
759
760
761
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
762
763
			case 'e' :
			element = strdup(optarg);
764
765
			break;

Thomas White's avatar
Thomas White committed
766
767
			case 2 :
			speaks = strdup(optarg);
768
769
			break;

Thomas White's avatar
Thomas White committed
770
771
			case 3 :
			scellr = strdup(optarg);
772
773
			break;

Thomas White's avatar
Thomas White committed
774
			case 4 :
775
776
777
			min_gradient = strtof(optarg, NULL);
			break;

Thomas White's avatar
Thomas White committed
778
			case 5 :
Thomas White's avatar
Thomas White committed
779
780
781
782
			stream_flags = parse_stream_flags(optarg);
			if ( stream_flags < 0 ) return 1;
			break;

Thomas White's avatar
Thomas White committed
783
			case 6 :
784
785
786
787
788
789
790
791
			cpu_num = strtol(optarg, &endptr, 10);
			if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
				ERROR("Invalid number of CPUs ('%s')\n",
				      optarg);
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
792
			case 7 :
793
794
795
796
797
798
799
800
801
802
803
804
805
			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;

Thomas White's avatar
Thomas White committed
806
			case 8 :
807
808
809
810
811
812
813
814
815
816
817
818
			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;

Thomas White's avatar
Thomas White committed
819
			case 9 :
820
821
822
			hdf5_peak_path = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
823
			case 10 :
824
825
826
			add_copy_hdf5_field(copyme, optarg);
			break;

Thomas White's avatar
Thomas White committed
827
828
829
830
831
832
			case 11 :
			min_snr = strtof(optarg, NULL);
			break;

			case 12 :
			min_int_snr = strtof(optarg, NULL);
Thomas White's avatar
Thomas White committed
833
834
			break;

Thomas White's avatar
Thomas White committed
835
836
837
838
			case 13 :
			toler = strdup(optarg);
			break;

839
840
841
842
			case 14 :
			intrad = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
843
844
845
846
			case 0 :
			break;

			default :
Thomas White's avatar
Thomas White committed
847
			return 1;
Thomas White's avatar
Thomas White committed
848

Thomas White's avatar
Thomas White committed
849
850
851
852
		}

	}

853
854
855
856
857
858
	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
859
860
861
862
863
864
865
866
867
	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
868
		ERROR("Failed to open input file '%s'\n", filename);
Thomas White's avatar
Thomas White committed
869
870
		return 1;
	}
Thomas White's avatar
Thomas White committed
871
	free(filename);
Thomas White's avatar
Thomas White committed
872

873
874
875
876
877
878
879
880
	if ( outfile == NULL ) {
		outfile = strdup("-");
	}
	if ( strcmp(outfile, "-") == 0 ) {
		ofh = stdout;
	} else {
		ofh = fopen(outfile, "w");
	}
Chuck's avatar
Chuck committed
881
//printf("***************%s\n",outfile);
882
883
884
885
886
	if ( ofh == NULL ) {
		ERROR("Failed to open output file '%s'\n", outfile);
		return 1;
	}

887
888
889
890
	if ( hdf5_peak_path == NULL ) {
		hdf5_peak_path = strdup("/processing/hitfinder/peakinfo");
	}

891
892
893
894
895
896
897
898
899
900
901
902
903
	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
904
	free(speaks);
905

Thomas White's avatar
Thomas White committed
906
907
908
909
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

910
	if ( prefix == NULL ) {
Thomas White's avatar
Thomas White committed
911
		prefix = strdup("");
912
	} else {
913
914
915
		if ( config_checkprefix ) {
			prefix = check_prefix(prefix);
		}
916
917
	}

Chuck's avatar
Chuck committed
918
	if ( nProcesses == 0 ) {
919
920
921
922
		ERROR("Invalid number of threads.\n");
		return 1;
	}

923
924
	if ( (indm_str == NULL) ||
	     ((indm_str != NULL) && (strcmp(indm_str, "none") == 0)) ) {
925
		STATUS("Not indexing anything.\n");
926
		indexer_needs_cell = 0;
927
928
929
		reduction_needs_cell = 0;
		indm = NULL;
		cellr = CELLR_NONE;
Thomas White's avatar
Thomas White committed
930
	} else {
931
932
933
934
935
936
937
938
		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
939
940
			indm = build_indexer_list(indm_str,
			                          &indexer_needs_cell);
941
942
943
944
945
			if ( indm == NULL ) {
				ERROR("Invalid indexer list '%s'\n", indm_str);
				return 1;
			}
			free(indm_str);
Thomas White's avatar
Thomas White committed
946
		}
Thomas White's avatar
Thomas White committed
947

948
949
950
951
952
953
954
		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 {
955
956
957
958
959
960
961
962
963
			int err;
			cellr = parse_cell_reduction(scellr, &err,
			                             &reduction_needs_cell);
			if ( err ) {
				ERROR("Unrecognised cell reduction '%s'\n",
			              scellr);
				return 1;
			}
			free(scellr);
964
		}
965
966
	}

967
968
969
	/* No indexing -> no reduction */
	if ( indm == NULL ) reduction_needs_cell = 0;

970
971
	if ( toler != NULL ) {
		int ttt;
Thomas White's avatar
Thomas White committed
972
973
		ttt = sscanf(toler, "%f,%f,%f,%f",
		             &tols[0], &tols[1], &tols[2], &tols[3] );
974
975
976
977
978
979
		if ( ttt != 4 ) {
			ERROR("Invalid parameters for '--tolerance'\n");
			return 1;
		}
	}

980
981
982
983
984
985
986
987
988
989
990
991
992
	if ( intrad != NULL ) {
		int r;
		r = sscanf(intrad, "%f,%f,%f", &ir_inn, &ir_mid, &ir_out);
		if ( r != 3 ) {
			ERROR("Invalid parameters for '--int-radius'\n");
			return 1;
		}
	} else {
		STATUS("WARNING: You did not specify --int-radius.\n");
		STATUS("WARNING: I will use the default values, which are"
		       " probably not appropriate for your patterns.\n");
	}

Thomas White's avatar
Thomas White committed
993
994
995
996
997
998
999
1000
	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);
For faster browsing, not all history is shown. View entire blame