/* * indexamajig.c * * Index patterns, output hkl+intensity etc. * * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: * 2010-2017 Thomas White * 2011 Richard Kirian * 2012 Lorenzo Galli * 2012 Chunhong Yoon * 2017 Valerio Mariani * 2017-2018 Yaroslav Gevorkov * * 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 . * */ #ifdef HAVE_CONFIG_H #include #endif #include #include #include #include #include #include #include #include #include #include #include #include "utils.h" #include "hdf5-file.h" #include "index.h" #include "peaks.h" #include "detector.h" #include "filters.h" #include "thread-pool.h" #include "geometry.h" #include "stream.h" #include "reflist-utils.h" #include "cell-utils.h" #include "integration.h" #include "taketwo.h" #include "im-sandbox.h" #include "image.h" static void show_help(const char *s) { printf("Syntax: %s [options]\n\n", s); printf( "Index and integrate snapshot diffraction images.\n\n" " -h, --help Display this help message\n" " --version Print CrystFEL version number and exit\n" "\nBasic options:\n\n" " -i, --input= List of images to process.\n" " -o, --output= Output stream filename\n" " -g --geometry= Detector geometry filename\n" " --basename Remove the directory parts of the filenames\n" " -x, --prefix=

Prefix filenames from input file with

\n" " --no-check-prefix Don't attempt to correct the --prefix\n" " -j Run analyses in parallel Default 1\n" " --highres= Absolute resolution cutoff in Angstroms\n" " --profile Show timing data for performance monitoring\n" " --temp-dir= Put the temporary folder under \n" " --wait-for-file= Time to wait for each file before processing\n" " --zmq-msgpack Receive data in MessagePack format over ZMQ\n" " --no-image-data Do not load image data (from ZMQ)\n" "\nPeak search options:\n\n" " --peaks= Peak search method. Default: zaef\n" " (zaef,peakfinder8,peakfinder9,hdf5,cxi,msgpack)\n" " --peak-radius= Integration radii for peak search\n" " --min-peaks= Minimum number of peaks for indexing\n" " --hdf5-peaks=

Find peaks table in HDF5 file here\n" " Default: /processing/hitfinder/peakinfo\n" " --median-filter= Apply a median filter to the image data\n" " Default: 0 (no filter)\n" " --filter-noise Apply noise filter to image data\n" " --threshold= Threshold for peak detection\n" " (zaef,peakfinder8 only) Default: 800\n" " --min-squared-gradient=\n" " Minimum squared gradient\n" " (zaef only) Default: 100,000\n" " --min-snr= Minimum signal/noise ratio for peaks\n" " (zaef,peakfinder8, peakfinder9 only) Default: 5\n" " --min-pix-count= Minimum number of pixels per peak\n" " (peakfinder8 only) Default: 2\n" " --max-pix-count= Maximum number of pixels per peak\n" " (peakfinder8 only) Default: 200\n" " --local-bg-radius= Radius (pixels) for local background estimation\n" " (peakfinder8, peakfinder9 only) Default: 3\n" " --min-res= Minimum resolution for peak search (in pixels)\n" " (peakfinder8 only) Default: 0\n" " --max-res= Maximum resolution for peak search (in pixels)\n" " (peakfinder8 only) Default: 1200\n" " --min-snr-biggest-pix=\n" " Minimum snr of the biggest pixel in the peak\n" " (peakfinder9 only)\n" " --min-snr-peak-pix=\n" " Minimum snr of a peak pixel (peakfinder9 only)\n" " --min-sig= Minimum standard deviation of the background\n" " (peakfinder9 only)\n" " --min-peak-over-neighbour=\n" " Just for speed. Biggest pixel in peak must be n\n" " higher than this (peakfinder9 only).\n" " --no-use-saturated Reject saturated peaks\n" " --no-revalidate Don't re-integrate and check HDF5 peaks\n" " --no-half-pixel-shift\n" " Don't offset the HDF5 peak locations by 0.5 px\n" " --check-hdf5-snr Check SNR for peaks from hdf5 or cxi (see --min-snr)\n" "\nIndexing options:\n\n" " --indexing= Indexing method list, comma separated\n" " -p, --pdb= Unit cell file (PDB or CrystFEL unit cell format)\n" " Default: 'molecule.pdb'\n" " --tolerance= Tolerances for cell comparison\n" " Default: 5,5,5,1.5,1.5,1.5\n" " --no-check-cell Don't check lattice parameters against input cell\n" " --no-cell-combinations\n" " Don't use axis combinations when checking cell\n" " --multi Repeat indexing to index multiple hits\n" " --no-retry Don't repeat indexing to increase indexing rate\n" " --no-refine Skip the prediction refinement step\n" " --no-check-peaks Don't check that most of the peaks can be accounted\n" " for by the indexing solution\n" "\n" " --taketwo-member-threshold\n" " Minimum number of members in network\n" " --taketwo-len-tolerance\n" " Reciprocal space length tolerance (1/A)\n" " --taketwo-angle-tolerance\n" " Reciprocal space angle tolerance (in degrees)\n" " --taketwo-trace-tolerance\n" " Rotation matrix equivalence tolerance (in degrees)\n" "\n" " --felix-domega Degree range of omega (moscaicity) to consider.\n" " Default: 2\n" " --felix-fraction-max-visits\n" " Cutoff for minimum fraction of the max visits.\n" " Default: 0.75\n" " --felix-max-internal-angle\n" " Cutoff for maximum internal angle between observed\n" " spots and predicted spots. Default: 0.25\n" " --felix-max-uniqueness\n" " Cutoff for maximum fraction of found spots which\n" " can belong to other crystallites. Default: 0.5\n" " --felix-min-completeness\n" " Cutoff for minimum fraction of projected spots\n" " found in the pattern. Default: 0.001\n" " --felix-min-visits\n" " Cutoff for minimum number of voxel visits.\n" " Default: 15\n" " --felix-num-voxels Number of voxels for Rodrigues space search\n" " Default: 100\n" " --felix-sigma The sigma of the 2theta, eta and omega angles.\n" " Default: 0.2\n" " --felix-tthrange-max Maximum 2theta to consider for indexing (degrees)\n" " Default: 30\n" " --felix-tthrange-min Minimum 2theta to consider for indexing (degrees)\n" " Default: 0\n" "\n" " --xgandalf-sampling-pitch\n" " Sampling pitch: 0 (loosest) to 4 (most dense)\n" " or with secondary Miller indices: 5 (loosest) to\n" " 7 (most dense). Default: 6\n" " --xgandalf-grad-desc-iterations\n" " Gradient descent iterations: 0 (few) to 5 (many)\n" " Default: 4\n" " --xgandalf-fast-execution Shortcut to set\n" " --xgandalf-sampling-pitch=2\n" " --xgandalf-grad-desc-iterations=3\n" " --xgandalf-tolerance Relative tolerance of the lattice vectors.\n" " Default is 0.02\n" " --xgandalf-no-deviation-from-provided-cell\n" " Force the fitted cell to have the same lattice\n" " parameters as the provided one\n" " --xgandalf-min-lattice-vector-length\n" " Minimum possible lattice vector length in A.\n" " Default: 30 A\n" " --xgandalf-max-lattice-vector-length\n" " Maximum possible lattice vector length in A.\n" " Default: 250 A\n" " --xgandalf-max-peaks\n" " Maximum number of peaks used for indexing.\n" " All peaks are used for refinement.\n" " Default: 250\n" "\n" "\nIntegration options:\n\n" " --integration= Integration method (rings,prof2d)-(cen,nocen)\n" " Default: rings-nocen\n" " --fix-profile-radius Fix the reciprocal space profile radius for spot\n" " prediction (default: automatically determine\n" " --fix-bandwidth Set the bandwidth for spot prediction\n" " --fix-divergence Set the divergence (full angle) for spot prediction\n" " --int-radius= Set the integration radii. Default: 4,5,7.\n" " --int-diag= Show debugging information about reflections\n" " --push-res= Integrate higher than apparent resolution cutoff\n" " --overpredict Over-predict reflections (for post-refinement)\n" "\nOutput options:\n\n" " --no-non-hits-in-stream\n" " Do not include non-hit frames in the stream\n" " (see --min-peaks)\n" " --copy-hdf5-field= Copy the value of HDF5 field into the stream\n" " --no-peaks-in-stream Do not record peak search results in the stream\n" " --no-refls-in-stream Do not record integrated reflections in the stream\n" " --serial-start Start the serial numbers in the stream here\n" "\nHistorical options:\n\n" " --no-sat-corr Don't correct values of saturated peaks\n" ); } static void add_geom_beam_stuff_to_field_list(struct imagefile_field_list *copyme, struct detector *det, struct beam_params *beam) { int i; for ( i=0; in_panels; i++ ) { struct panel *p = &det->panels[i]; if ( p->clen_from != NULL ) { add_imagefile_field(copyme, p->clen_from); } } if ( beam->photon_energy_from != NULL ) { add_imagefile_field(copyme, beam->photon_energy_from); } } int main(int argc, char *argv[]) { int c; char *filename = NULL; char *outfile = NULL; FILE *fh; Stream *st; int config_checkprefix = 1; int config_basename = 0; int integrate_saturated = 0; char *indm_str = NULL; char *cellfile = NULL; char *prefix = NULL; char *speaks = NULL; char *toler = NULL; int n_proc = 1; struct index_args iargs; char *intrad = NULL; char *pkrad = NULL; char *int_str = NULL; char *temp_location = NULL; /* e.g. /tmp */ char *tmpdir; /* e.g. /tmp/indexamajig.12345 */ char *rn; /* e.g. /home/taw/indexing */ int r; char *int_diag = NULL; char *geom_filename = NULL; struct beam_params beam; int have_push_res = 0; char *command_line_peak_path = NULL; int if_refine = 1; int if_nocomb = 0; int if_nocheck = 0; int if_peaks = 1; int if_multi = 0; int if_retry = 1; int serial_start = 1; char *spectrum_fn = NULL; int zmq = 0; char *zmq_address = NULL; /* Defaults */ iargs.cell = NULL; iargs.noisefilter = 0; iargs.median_filter = 0; iargs.satcorr = 1; iargs.tols[0] = 0.05; iargs.tols[1] = 0.05; iargs.tols[2] = 0.05; iargs.tols[3] = 1.5; iargs.tols[4] = 1.5; iargs.tols[5] = 1.5; iargs.threshold = 800.0; iargs.min_sq_gradient = 100000.0; iargs.min_snr = 5.0; iargs.min_pix_count = 2; iargs.max_pix_count = 200; iargs.min_res = 0; iargs.max_res = 1200; iargs.local_bg_radius = 3; iargs.min_snr_biggest_pix = 7.0; /* peak finder 9 */ iargs.min_snr_peak_pix = 6.0; iargs.min_sig = 11.0; iargs.min_peak_over_neighbour = -INFINITY; iargs.check_hdf5_snr = 0; iargs.det = NULL; iargs.peaks = PEAK_ZAEF; iargs.beam = &beam; iargs.hdf5_peak_path = NULL; iargs.half_pixel_shift = 1; iargs.copyme = NULL; iargs.pk_inn = -1.0; iargs.pk_mid = -1.0; iargs.pk_out = -1.0; iargs.ir_inn = 4.0; iargs.ir_mid = 5.0; iargs.ir_out = 7.0; iargs.use_saturated = 1; iargs.no_revalidate = 0; iargs.stream_peaks = 1; iargs.stream_refls = 1; iargs.stream_nonhits = 1; iargs.int_diag = INTDIAG_NONE; iargs.copyme = new_imagefile_field_list(); iargs.min_peaks = 0; iargs.overpredict = 0; iargs.wait_for_file = 0; if ( iargs.copyme == NULL ) { ERROR("Couldn't allocate HDF5 field list.\n"); return 1; } iargs.ipriv = NULL; /* No default */ iargs.int_meth = integration_method("rings-nocen-nosat-nograd", NULL); iargs.push_res = 0.0; iargs.highres = +INFINITY; iargs.fix_profile_r = -1.0; iargs.fix_bandwidth = -1.0; iargs.fix_divergence = -1.0; iargs.profile = 0; iargs.no_image_data = 0; iargs.taketwo_opts.member_thresh = -1; iargs.taketwo_opts.len_tol = -1.0; iargs.taketwo_opts.angle_tol = -1.0; iargs.taketwo_opts.trace_tol = -1.0; iargs.xgandalf_opts.sampling_pitch = 6; iargs.xgandalf_opts.grad_desc_iterations = 4; iargs.xgandalf_opts.tolerance = 0.02; iargs.xgandalf_opts.no_deviation_from_provided_cell = 0; iargs.xgandalf_opts.minLatticeVectorLength_A = 30; iargs.xgandalf_opts.maxLatticeVectorLength_A = 250; iargs.xgandalf_opts.maxPeaksForIndexing = 250; iargs.felix_opts.ttmin = -1.0; iargs.felix_opts.ttmax = -1.0; iargs.felix_opts.min_visits = 0; iargs.felix_opts.min_completeness = -1.0; iargs.felix_opts.max_uniqueness = -1.0; iargs.felix_opts.n_voxels = 0; iargs.felix_opts.fraction_max_visits = -1.0; iargs.felix_opts.sigma = -1.0; iargs.felix_opts.domega = -1.0; iargs.felix_opts.max_internal_angle = -1.0; /* Long options */ const struct option longopts[] = { /* Options with long and short versions */ {"help", 0, NULL, 'h'}, {"version", 0, NULL, 'v'}, {"input", 1, NULL, 'i'}, {"output", 1, NULL, 'o'}, {"indexing", 1, NULL, 'z'}, {"geometry", 1, NULL, 'g'}, {"pdb", 1, NULL, 'p'}, {"prefix", 1, NULL, 'x'}, {"threshold", 1, NULL, 't'}, {"beam", 1, NULL, 'b'}, /* Long-only options with no arguments */ {"filter-noise", 0, &iargs.noisefilter, 1}, {"no-check-prefix", 0, &config_checkprefix, 0}, {"basename", 0, &config_basename, 1}, {"no-peaks-in-stream", 0, &iargs.stream_peaks, 0}, {"no-refls-in-stream", 0, &iargs.stream_refls, 0}, {"no-non-hits-in-stream", 0, &iargs.stream_nonhits, 0}, {"integrate-saturated",0, &integrate_saturated, 1}, {"no-use-saturated", 0, &iargs.use_saturated, 0}, {"no-revalidate", 0, &iargs.no_revalidate, 1}, {"check-hdf5-snr", 0, &iargs.check_hdf5_snr, 1}, {"profile", 0, &iargs.profile, 1}, {"no-half-pixel-shift",0, &iargs.half_pixel_shift, 0}, {"no-refine", 0, &if_refine, 0}, {"no-cell-combinations",0,&if_nocomb, 1}, {"no-check-cell", 0, &if_nocheck, 1}, {"no-cell-check", 0, &if_nocheck, 1}, {"check-peaks", 0, &if_peaks, 1}, {"no-check-peaks", 0, &if_peaks, 0}, {"no-retry", 0, &if_retry, 0}, {"retry", 0, &if_retry, 1}, {"no-multi", 0, &if_multi, 0}, {"multi", 0, &if_multi, 1}, {"overpredict", 0, &iargs.overpredict, 1}, {"zmq-msgpack", 0, &zmq, 1}, {"no-image-data", 0, &iargs.no_image_data, 1}, /* Long-only options which don't actually do anything */ {"no-sat-corr", 0, &iargs.satcorr, 0}, {"sat-corr", 0, &iargs.satcorr, 1}, {"no-check-hdf5-snr", 0, &iargs.check_hdf5_snr, 0}, {"use-saturated", 0, &iargs.use_saturated, 1}, /* Long-only options with arguments */ {"peaks", 1, NULL, 302}, {"cell-reduction", 1, NULL, 303}, {"min-gradient", 1, NULL, 304}, {"record", 1, NULL, 305}, {"cpus", 1, NULL, 306}, {"cpugroup", 1, NULL, 307}, {"cpuoffset", 1, NULL, 308}, {"hdf5-peaks", 1, NULL, 309}, {"copy-hdf5-field", 1, NULL, 310}, {"min-snr", 1, NULL, 311}, {"tolerance", 1, NULL, 313}, {"int-radius", 1, NULL, 314}, {"median-filter", 1, NULL, 315}, {"integration", 1, NULL, 316}, {"temp-dir", 1, NULL, 317}, {"int-diag", 1, NULL, 318}, {"push-res", 1, NULL, 319}, {"res-push", 1, NULL, 319}, /* compat */ {"peak-radius", 1, NULL, 320}, {"highres", 1, NULL, 321}, {"fix-profile-radius", 1, NULL, 322}, {"fix-bandwidth", 1, NULL, 323}, {"fix-divergence", 1, NULL, 324}, {"felix-options", 1, NULL, 325}, {"min-pix-count", 1, NULL, 326}, {"max-pix-count", 1, NULL, 327}, {"local-bg-radius", 1, NULL, 328}, {"min-res", 1, NULL, 329}, {"max-res", 1, NULL, 330}, {"min-peaks", 1, NULL, 331}, {"taketwo-member-threshold", 1, NULL, 332}, {"taketwo-member-thresh", 1, NULL, 332}, /* compat */ {"taketwo-len-tolerance", 1, NULL, 333}, {"taketwo-len-tol", 1, NULL, 333}, /* compat */ {"taketwo-angle-tolerance", 1, NULL, 334}, {"taketwo-angle-tol", 1, NULL, 334}, /* compat */ {"taketwo-trace-tolerance", 1, NULL, 335}, {"taketwo-trace-tol", 1, NULL, 335}, /* compat */ {"felix-tthrange-min", 1, NULL, 336}, {"felix-tthrange-max", 1, NULL, 337}, {"felix-min-visits", 1, NULL, 338}, {"felix-min-completeness", 1, NULL, 339}, {"felix-max-uniqueness", 1, NULL, 340}, {"felix-num-voxels", 1, NULL, 341}, {"felix-fraction-max-visits",1, NULL, 342}, {"felix-sigma", 1, NULL, 343}, {"serial-start", 1, NULL, 344}, {"felix-domega", 1, NULL, 345}, {"felix-max-internal-angle", 1, NULL, 346}, {"min-snr-biggest-pix", 1, NULL, 347}, {"min-snr-peak-pix", 1, NULL, 348}, {"min-sig", 1, NULL, 349}, {"min-peak-over-neighbour", 1, NULL, 350}, {"xgandalf-sampling-pitch", 1, NULL, 351}, {"xgandalf-sps", 1, NULL, 351}, {"xgandalf-grad-desc-iterations", 1, NULL, 352}, {"xgandalf-gdis", 1, NULL, 352}, {"xgandalf-tolerance", 1, NULL, 353}, {"xgandalf-tol", 1, NULL, 353}, {"xgandalf-no-deviation-from-provided-cell", 0, NULL, 354}, {"xgandalf-ndfpc", 0, NULL, 354}, {"xgandalf-min-lattice-vector-length", 1, NULL, 355}, {"xgandalf-min-lvl", 1, NULL, 355}, {"xgandalf-max-lattice-vector-length", 1, NULL, 356}, {"xgandalf-max-lvl", 1, NULL, 356}, {"spectrum-file", 1, NULL, 357}, {"wait-for-file", 1, NULL, 358}, {"min-squared-gradient",1,NULL, 359}, {"min-sq-gradient", 1, NULL, 359}, /* compat */ {"xgandalf-fast-execution", 0, NULL, 360}, {"xgandalf-max-peaks", 1, NULL, 361}, {0, 0, NULL, 0} }; /* Short options */ while ((c = getopt_long(argc, argv, "hi:o:z:p:x:j:g:t:vb:", longopts, NULL)) != -1) { switch (c) { case 'h' : show_help(argv[0]); return 0; case 'v' : printf("CrystFEL: " CRYSTFEL_VERSIONSTRING "\n"); printf(CRYSTFEL_BOILERPLATE"\n"); return 0; case 'b' : ERROR("WARNING: This version of CrystFEL no longer " "uses beam files. Please remove the beam file " "from your indexamajig command line.\n"); return 1; case 'i' : filename = strdup(optarg); break; case 'o' : outfile = strdup(optarg); break; case 'z' : indm_str = strdup(optarg); break; case 'p' : cellfile = strdup(optarg); break; case 'x' : prefix = strdup(optarg); break; case 'j' : n_proc = atoi(optarg); break; case 'g' : geom_filename = optarg; break; case 't' : iargs.threshold = strtof(optarg, NULL); break; case 302 : speaks = strdup(optarg); break; case 303 : ERROR("The option '--cell-reduction' is no longer " "used.\n" "The complete indexing behaviour is now " "controlled using '--indexing'.\n" "See 'man indexamajig' for details of the " "available methods.\n"); return 1; case 304 : iargs.min_sq_gradient = strtof(optarg, NULL); ERROR("Recommend using --min-squared-gradient instead " "of --min-gradient.\n"); break; case 305 : ERROR("The option '--record' is no longer used.\n" "Use '--no-peaks-in-stream' and" "'--no-refls-in-stream' if you need to control" "the contents of the stream.\n"); return 1; case 306 : case 307 : case 308 : ERROR("The options --cpus, --cpugroup and --cpuoffset" " are no longer used by indexamajig.\n"); break; case 309 : free(command_line_peak_path); command_line_peak_path = strdup(optarg); break; case 310 : add_imagefile_field(iargs.copyme, optarg); break; case 311 : iargs.min_snr = strtof(optarg, NULL); break; case 313 : toler = strdup(optarg); break; case 314 : intrad = strdup(optarg); break; case 315 : iargs.median_filter = atoi(optarg); break; case 316 : int_str = strdup(optarg); break; case 317 : temp_location = strdup(optarg); break; case 318 : int_diag = strdup(optarg); break; case 319 : if ( sscanf(optarg, "%f", &iargs.push_res) != 1 ) { ERROR("Invalid value for --push-res\n"); return 1; } iargs.push_res *= 1e9; /* nm^-1 -> m^-1 */ have_push_res = 1; break; case 320 : pkrad = strdup(optarg); break; case 321 : if ( sscanf(optarg, "%f", &iargs.highres) != 1 ) { ERROR("Invalid value for --highres\n"); return 1; } /* A -> m^-1 */ iargs.highres = 1.0 / (iargs.highres/1e10); break; case 322 : if ( sscanf(optarg, "%f", &iargs.fix_profile_r) != 1 ) { ERROR("Invalid value for " "--fix-profile-radius\n"); return 1; } break; case 323 : if ( sscanf(optarg, "%f", &iargs.fix_bandwidth) != 1 ) { ERROR("Invalid value for --fix-bandwidth\n"); return 1; } break; case 324 : if ( sscanf(optarg, "%f", &iargs.fix_divergence) != 1 ) { ERROR("Invalid value for --fix-divergence\n"); return 1; } break; case 325 : ERROR("--felix-options is no longer used.\n"); ERROR("See --help for the new Felix options.\n"); return 1; case 326: iargs.min_pix_count = atoi(optarg); break; case 327: iargs.max_pix_count = atoi(optarg); break; case 328: iargs.local_bg_radius = atoi(optarg); break; case 329: iargs.min_res = atoi(optarg); break; case 330: iargs.max_res = atoi(optarg); break; case 331: iargs.min_peaks = atoi(optarg); break; case 332: if ( sscanf(optarg, "%i", &iargs.taketwo_opts.member_thresh) != 1 ) { ERROR("Invalid value for --taketwo-member-threshold\n"); return 1; } break; case 333: if ( sscanf(optarg, "%lf", &iargs.taketwo_opts.len_tol) != 1 ) { ERROR("Invalid value for --taketwo-len-tolerance\n"); return 1; } /* Convert to m^-1 */ iargs.taketwo_opts.len_tol *= 1e10; break; case 334: if ( sscanf(optarg, "%lf", &iargs.taketwo_opts.angle_tol) != 1 ) { ERROR("Invalid value for --taketwo-angle-tolerance\n"); return 1; } /* Convert to radians */ iargs.taketwo_opts.angle_tol = deg2rad(iargs.taketwo_opts.angle_tol); break; case 335: if ( sscanf(optarg, "%lf", &iargs.taketwo_opts.trace_tol) != 1 ) { ERROR("Invalid value for --taketwo-trace-tolerance\n"); return 1; } /* Convert to radians */ iargs.taketwo_opts.trace_tol = deg2rad(iargs.taketwo_opts.trace_tol); break; case 336: if ( sscanf(optarg, "%lf", &iargs.felix_opts.ttmin) != 1 ) { ERROR("Invalid value for --felix-tthrange-min\n"); return 1; } iargs.felix_opts.ttmin = deg2rad(iargs.felix_opts.ttmin); break; case 337: if ( sscanf(optarg, "%lf", &iargs.felix_opts.ttmax) != 1 ) { ERROR("Invalid value for --felix-tthrange-max\n"); return 1; } iargs.felix_opts.ttmax = deg2rad(iargs.felix_opts.ttmax); break; case 338: if ( sscanf(optarg, "%i", &iargs.felix_opts.min_visits) != 1 ) { ERROR("Invalid value for --felix-min-visits\n"); return 1; } break; case 339: if ( sscanf(optarg, "%lf", &iargs.felix_opts.min_completeness) != 1 ) { ERROR("Invalid value for --felix-min-completeness\n"); return 1; } break; case 340: if ( sscanf(optarg, "%lf", &iargs.felix_opts.max_uniqueness) != 1 ) { ERROR("Invalid value for --felix-max-uniqueness\n"); return 1; } break; case 341: if ( sscanf(optarg, "%i", &iargs.felix_opts.n_voxels) != 1 ) { ERROR("Invalid value for --felix-num-voxels\n"); return 1; } break; case 342: if ( sscanf(optarg, "%lf", &iargs.felix_opts.fraction_max_visits) != 1 ) { ERROR("Invalid value for --felix-fraction-max-visits\n"); return 1; } break; case 343: if ( sscanf(optarg, "%lf", &iargs.felix_opts.sigma) != 1 ) { ERROR("Invalid value for --felix-sigma\n"); return 1; } break; case 344: if ( sscanf(optarg, "%i", &serial_start) != 1 ) { ERROR("Invalid value for --serial-start\n"); return 1; } break; case 345: if ( sscanf(optarg, "%lf", &iargs.felix_opts.domega) != 1 ) { ERROR("Invalid value for --felix-domega\n"); return 1; } break; case 346: if ( sscanf(optarg, "%lf", &iargs.felix_opts.max_internal_angle) != 1 ) { ERROR("Invalid value for --felix-max-internal-angle\n"); return 1; } break; case 347: iargs.min_snr_biggest_pix = strtof(optarg, NULL); break; case 348: iargs.min_snr_peak_pix = strtof(optarg, NULL); break; case 349: iargs.min_sig = strtof(optarg, NULL); break; case 350: iargs.min_peak_over_neighbour = strtof(optarg, NULL); break; case 351: if (sscanf(optarg, "%u", &iargs.xgandalf_opts.sampling_pitch) != 1) { ERROR("Invalid value for --xgandalf-sampling-pitch\n"); return 1; } break; case 352: if (sscanf(optarg, "%u", &iargs.xgandalf_opts.grad_desc_iterations) != 1) { ERROR("Invalid value for --xgandalf-grad-desc-iterations\n"); return 1; } break; case 353: if (sscanf(optarg, "%f", &iargs.xgandalf_opts.tolerance) != 1) { ERROR("Invalid value for --xgandalf-tolerance\n"); return 1; } break; case 354: iargs.xgandalf_opts.no_deviation_from_provided_cell = 1; break; case 355: if (sscanf(optarg, "%f", &iargs.xgandalf_opts.minLatticeVectorLength_A) != 1) { ERROR("Invalid value for " "--xgandalf-min-lattice-vector-length\n"); return 1; } break; case 356: if (sscanf(optarg, "%f", &iargs.xgandalf_opts.maxLatticeVectorLength_A) != 1) { ERROR("Invalid value for " "--xgandalf-max-lattice-vector-length\n"); return 1; } break; case 357: spectrum_fn = strdup(optarg); break; case 358: if (sscanf(optarg, "%d", &iargs.wait_for_file) != 1) { ERROR("Invalid value for --wait-for-file\n"); return 1; } break; case 359 : iargs.min_sq_gradient = strtof(optarg, NULL); break; case 360: iargs.xgandalf_opts.sampling_pitch = 2; iargs.xgandalf_opts.grad_desc_iterations = 3; break; case 361: if (sscanf(optarg, "%i", &iargs.xgandalf_opts.maxPeaksForIndexing) != 1) { ERROR("Invalid value for --xgandalf-max-peaks\n"); return 1; } break; case 0 : break; case '?' : break; default : ERROR("Unhandled option '%c'\n", c); break; } } /* Check for minimal information */ if ( filename == NULL ) { ERROR("You need to provide the input filename (use -i)\n"); return 1; } if ( geom_filename == NULL ) { ERROR("You need to specify the geometry filename (use -g)\n"); return 1; } if ( outfile == NULL ) { ERROR("You need to specify the output filename (use -o)\n"); return 1; } if ( temp_location == NULL ) { temp_location = strdup("."); } /* Open input */ if ( strcmp(filename, "-") == 0 ) { fh = stdin; } else { fh = fopen(filename, "r"); } if ( fh == NULL ) { ERROR("Failed to open input file '%s'\n", filename); return 1; } free(filename); /* Parse peak detection method */ 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 ) { iargs.peaks = PEAK_ZAEF; } else if ( strcmp(speaks, "peakfinder8") == 0 ) { iargs.peaks = PEAK_PEAKFINDER8; } else if ( strcmp(speaks, "hdf5") == 0 ) { iargs.peaks = PEAK_HDF5; } else if ( strcmp(speaks, "cxi") == 0 ) { iargs.peaks = PEAK_CXI; } else if ( strcmp(speaks, "peakfinder9") == 0 ) { iargs.peaks = PEAK_PEAKFINDER9; } else if ( strcmp(speaks, "msgpack") == 0 ) { iargs.peaks = PEAK_MSGPACK; } else { ERROR("Unrecognised peak detection method '%s'\n", speaks); return 1; } free(speaks); /* Check prefix (if given) */ if ( prefix == NULL ) { prefix = strdup(""); } else { if ( config_checkprefix ) { prefix = check_prefix(prefix); } } /* Check number of processes */ if ( n_proc == 0 ) { ERROR("Invalid number of processes.\n"); return 1; } /* Load detector geometry */ iargs.det = get_detector_geometry_2(geom_filename, iargs.beam, &iargs.hdf5_peak_path); if ( iargs.det == NULL ) { ERROR("Failed to read detector geometry from '%s'\n", geom_filename); return 1; } add_geom_beam_stuff_to_field_list(iargs.copyme, iargs.det, iargs.beam); /* If no peak path from geometry file, use these (but see later) */ if ( iargs.hdf5_peak_path == NULL ) { if ( iargs.peaks == PEAK_HDF5 ) { iargs.hdf5_peak_path = strdup("/processing/hitfinder/peakinfo"); } else if ( iargs.peaks == PEAK_CXI ) { iargs.hdf5_peak_path = strdup("/entry_1/result_1"); } } /* If an HDF5 peak path was given on the command line, use it */ if ( command_line_peak_path != NULL ) { free(iargs.hdf5_peak_path); iargs.hdf5_peak_path = command_line_peak_path; } /* Parse integration method */ if ( int_str != NULL ) { int err; iargs.int_meth = integration_method(int_str, &err); if ( err ) { ERROR("Invalid integration method '%s'\n", int_str); return 1; } free(int_str); } if ( integrate_saturated ) { /* Option provided for backwards compatibility */ iargs.int_meth |= INTEGRATION_SATURATED; } if ( have_push_res && !(iargs.int_meth & INTEGRATION_RESCUT) ) { ERROR("WARNING: You used --push-res, but not -rescut, " "therefore --push-res will have no effect.\n"); ERROR("WARNING: Add --integration=rings-rescut or " "--integration=prof2d-rescut.\n"); } /* Parse unit cell tolerance */ if ( toler != NULL ) { int ttt; ttt = sscanf(toler, "%f,%f,%f,%f,%f,%f", &iargs.tols[0], &iargs.tols[1], &iargs.tols[2], &iargs.tols[3], &iargs.tols[4], &iargs.tols[5]); if ( ttt != 6 ) { ttt = sscanf(toler, "%f,%f,%f,%f", &iargs.tols[0], &iargs.tols[1], &iargs.tols[2], &iargs.tols[3]); if ( ttt != 4 ) { ERROR("Invalid parameters for '--tolerance'\n"); return 1; } iargs.tols[4] = iargs.tols[3]; iargs.tols[5] = iargs.tols[3]; } free(toler); /* Percent to fraction */ iargs.tols[0] /= 100.0; iargs.tols[1] /= 100.0; iargs.tols[2] /= 100.0; iargs.tols[3] = deg2rad(iargs.tols[3]); iargs.tols[4] = deg2rad(iargs.tols[4]); iargs.tols[5] = deg2rad(iargs.tols[5]); } /* Parse integration radii */ if ( intrad != NULL ) { int r; r = sscanf(intrad, "%f,%f,%f", &iargs.ir_inn, &iargs.ir_mid, &iargs.ir_out); if ( r != 3 ) { ERROR("Invalid parameters for '--int-radius'\n"); return 1; } free(intrad); } 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"); } /* Parse peak radii (used for peak detection) */ if ( pkrad != NULL ) { int r; r = sscanf(pkrad, "%f,%f,%f", &iargs.pk_inn, &iargs.pk_mid, &iargs.pk_out); if ( r != 3 ) { ERROR("Invalid parameters for '--peak-radius'\n"); return 1; } free(pkrad); } if ( iargs.pk_inn < 0.0 ) { iargs.pk_inn = iargs.ir_inn; iargs.pk_mid = iargs.ir_mid; iargs.pk_out = iargs.ir_out; } /* Load unit cell (if given) */ if ( cellfile != NULL ) { iargs.cell = load_cell_from_file(cellfile); if ( iargs.cell == NULL ) { ERROR("Couldn't read unit cell (from %s)\n", cellfile); return 1; } free(cellfile); } else { iargs.cell = NULL; } /* Load spectrum from file if given */ if ( spectrum_fn != NULL ) { iargs.spectrum = spectrum_load(spectrum_fn); if ( iargs.spectrum == NULL ) { ERROR("Couldn't read spectrum (from %s)\n", spectrum_fn); return 1; } free(spectrum_fn); } else { iargs.spectrum = NULL; } /* Parse integration diagnostic */ if ( int_diag != NULL ) { int r; signed int h, k, l; if ( strcmp(int_diag, "random") == 0 ) { iargs.int_diag = INTDIAG_RANDOM; } if ( strcmp(int_diag, "all") == 0 ) { iargs.int_diag = INTDIAG_ALL; } if ( strcmp(int_diag, "negative") == 0 ) { iargs.int_diag = INTDIAG_NEGATIVE; } if ( strcmp(int_diag, "implausible") == 0 ) { iargs.int_diag = INTDIAG_IMPLAUSIBLE; } if ( strcmp(int_diag, "strong") == 0 ) { iargs.int_diag = INTDIAG_STRONG; } r = sscanf(int_diag, "%i,%i,%i", &h, &k, &l); if ( r == 3 ) { iargs.int_diag = INTDIAG_INDICES; iargs.int_diag_h = h; iargs.int_diag_k = k; iargs.int_diag_l = l; } if ( (iargs.int_diag == INTDIAG_NONE) && (strcmp(int_diag, "none") != 0) ) { ERROR("Invalid value for --int-diag.\n"); return 1; } free(int_diag); } tmpdir = create_tempdir(temp_location); if ( tmpdir == NULL ) return 1; /* Change into temporary folder, temporarily, to control the crap * dropped by indexing programs during setup */ rn = getcwd(NULL, 0); r = chdir(tmpdir); if ( r ) { ERROR("Failed to chdir to temporary folder: %s\n", strerror(errno)); return 1; } if ( indm_str == NULL ) { STATUS("No indexing methods specified. I will try to "); STATUS("automatically detect the available methods.\n"); STATUS("To disable auto-detection of indexing methods, specify "); STATUS("which methods to use with --indexing=.\n"); STATUS("Use --indexing=none to disable indexing and integration.\n"); indm_str = detect_indexing_methods(iargs.cell); } /* Prepare the indexing system */ if ( indm_str == NULL ) { ERROR("No indexing method specified, and no usable indexing "); ERROR("methods auto-detected.\n"); ERROR("Install some indexing programs (mosflm,dirax etc), or "); ERROR("try again with --indexing=none.\n"); return 1; } else if ( strcmp(indm_str, "none") == 0 ) { STATUS("Indexing/integration disabled.\n"); if ( iargs.cell != NULL ) { STATUS("Ignoring your unit cell.\n"); } iargs.ipriv = NULL; } else { IndexingFlags flags = 0; if ( iargs.cell != NULL ) { STATUS("This is what I understood your unit cell to be:\n"); cell_print(iargs.cell); } else { STATUS("No reference unit cell provided.\n"); } if ( if_nocomb ) { flags |= INDEXING_CHECK_CELL_AXES; } else { flags |= INDEXING_CHECK_CELL_COMBINATIONS; } if ( if_nocheck ) { flags &= ~INDEXING_CHECK_CELL_AXES; flags &= ~INDEXING_CHECK_CELL_COMBINATIONS; } if ( if_refine ) { flags |= INDEXING_REFINE; } if ( if_peaks ) { flags |= INDEXING_CHECK_PEAKS; } if ( if_multi ) { flags |= INDEXING_MULTI; } if ( if_retry ) { flags |= INDEXING_RETRY; } iargs.ipriv = setup_indexing(indm_str, iargs.cell, iargs.det, iargs.tols, flags, &iargs.taketwo_opts, &iargs.xgandalf_opts, &iargs.felix_opts); if ( iargs.ipriv == NULL ) { ERROR("Failed to set up indexing system\n"); return 1; } } /* Change back to where we were before. Sandbox code will create * worker subdirs inside the temporary folder, and process_image will * change into them. */ r = chdir(rn); if ( r ) { ERROR("Failed to chdir: %s\n", strerror(errno)); return 1; } free(rn); /* Open output stream */ st = open_stream_for_write_4(outfile, geom_filename, iargs.cell, argc, argv, indm_str); if ( st == NULL ) { ERROR("Failed to open stream '%s'\n", outfile); return 1; } free(outfile); free(indm_str); gsl_set_error_handler_off(); if ( zmq ) { char line[1024]; char *rval; rval = fgets(line, 1024, fh); if ( rval == NULL ) { ERROR("Failed to read ZMQ server/port from input.\n"); return 1; } chomp(line); zmq_address = strdup(line); /* In future, read multiple addresses and hand them out * evenly to workers */ } r = create_sandbox(&iargs, n_proc, prefix, config_basename, fh, st, tmpdir, serial_start, zmq_address); free_imagefile_field_list(iargs.copyme); cell_free(iargs.cell); free(iargs.beam->photon_energy_from); free(prefix); free(temp_location); free(tmpdir); free_detector_geometry(iargs.det); free(iargs.hdf5_peak_path); close_stream(st); cleanup_indexing(iargs.ipriv); if ( r ) { return 0; } else { return 1; /* No patterns processed */ } }