diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 10cc300c..0f57c359 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include "../reader/JFJochHDF5Reader.h" #include "../common/Logger.h" @@ -32,43 +33,87 @@ #include "../image_analysis/UpdateReflectionResolution.h" void print_usage(Logger &logger) { - logger.Info("Usage ./jfjoch_analysis {} "); - logger.Info("Options:"); - logger.Info(" -o Output file prefix (default: output)"); - logger.Info(" -N Number of threads (default: 1)"); - logger.Info(" -s Start image number (default: 0)"); - logger.Info(" -e End image number (default: all)"); - logger.Info(" -v Verbose output"); - logger.Info(""); + std::cout << "Usage ./jfjoch_analysis {} " << std::endl; + std::cout << "Options:" << std::endl; + std::cout << " -o, --output-prefix Output file prefix (default: output)" << std::endl; + std::cout << " -N, --threads Number of threads (default: 1)" << std::endl; + std::cout << " -s, --start-image Start image number (default: 0)" << std::endl; + std::cout << " -e, --end-image End image number (default: all)" << std::endl; + std::cout << " -t, --stride Image stride (default: 1)" << std::endl; + std::cout << " -v, --verbose Verbose output" << std::endl; + std::cout << std::endl; - logger.Info(" Spot finding"); - logger.Info(" -T Noise sigma level for spot finding (default: 3.0)"); - logger.Info(" -t Photon count threshold for spot finding (default: 10)"); - logger.Info(" -d High resolution limit for spot finding (default: 1.5)"); - logger.Info(" -c Max spot count (default: 250)"); - logger.Info(""); - logger.Info(" Indexing"); - logger.Info(" -R[num] Rotation indexing (optional: min angular range deg)"); - logger.Info(" -X Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)"); - logger.Info(" -F Use FFT indexing algorithm (shortcut for -XFFT)"); - logger.Info(" -S Space group number - used for both indexing and scaling"); - logger.Info( - " -C Fix reference unit cell: -C\"a,b,c,alpha,beta,gamma\" (comma-separated, no spaces; quotes optional)"); - logger.Info(" -x No least-square beam center refinement"); + std::cout << " Spot finding" << std::endl; + std::cout << " --spot-sigma Noise sigma level for spot finding (default: 3.0)" << std::endl; + std::cout << " --spot-threshold Photon count threshold for spot finding (default: 10)" << std::endl; + std::cout << " --spot-resolution High resolution limit for spot finding (default: 1.5)" << std::endl; + std::cout << " --max-spots Max spot count (default: 250)" << std::endl; + std::cout << std::endl; - logger.Info(""); - logger.Info(" Scaling and merging"); - logger.Info(" -D High resolution limit for scaling/merging (default: 0.0; no limit)"); - logger.Info(" -M Scale and merge (refine mosaicity) and write scaled.hkl + image.dat"); - logger.Info(" -P Partiality refinement fixed|rot|unity (default: fixed)"); - logger.Info(" -A Anomalous mode (don't merge Friedel pairs)"); - logger.Info(" -B Refine per image B-factor"); - logger.Info(" -w Refine image wedge during scaling with starting wedge value"); - logger.Info(" -p Minimum partiality to accept reflection (default: 0.02)"); - logger.Info(" -q Per-image CC limit in percent (default: no limit)"); - logger.Info(" -i Number of scaling iterations with no reference data (default: 3) "); + std::cout << " Indexing" << std::endl; + std::cout << " -R, --rotation-indexing[=num] Rotation indexing (optional: min angular range deg)" << std::endl; + std::cout << " -X, --indexing-algorithm Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl; + std::cout << " -F, --fft-indexing Use FFT indexing algorithm (shortcut for -X FFT)" << std::endl; + std::cout << " -S, --space-group Space group number - used for both indexing and scaling" << std::endl; + std::cout << " -C, --unit-cell Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl; + std::cout << " --no-beam-center-refine No least-square beam center refinement" << std::endl; + std::cout << std::endl; + + std::cout << " Scaling and merging" << std::endl; + std::cout << " --scale-merge Scale and merge (refine mosaicity) and write scaled.hkl + image.dat" << std::endl; + std::cout << " --partiality Partiality refinement fixed|rot|unity (default: fixed)" << std::endl; + std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl; + std::cout << " -B, --refine-bfactor Refine per image B-factor" << std::endl; + std::cout << " -w, --wedge[=num] Refine image wedge during scaling with starting wedge value" << std::endl; + std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl; + std::cout << " --min-image-cc Per-image CC limit in percent (default: no limit)" << std::endl; + std::cout << " --scaling-iterations Number of scaling iterations with no reference data (default: 3)" << std::endl; + std::cout << " -z, --reference-mtz Reference MTZ file" << std::endl; } +enum { + OPT_SPOT_SIGMA = 1000, + OPT_SPOT_THRESHOLD, + OPT_SPOT_RESOLUTION, + OPT_MAX_SPOTS, + OPT_NO_BEAM_CENTER_REFINE, + OPT_SCALE_MERGE, + OPT_MIN_PARTIALITY, + OPT_MIN_IMAGE_CC, + OPT_SCALING_ITERATIONS +}; + +static option long_options[] = { + {"verbose", no_argument, nullptr, 'v'}, + {"output-prefix", required_argument, nullptr, 'o'}, + {"threads", required_argument, nullptr, 'N'}, + {"start-image", required_argument, nullptr, 's'}, + {"end-image", required_argument, nullptr, 'e'}, +{"stride", required_argument, nullptr, 't'}, + {"rotation-indexing", optional_argument, nullptr, 'R'}, + {"indexing-algorithm", required_argument, nullptr, 'X'}, + {"unit-cell", required_argument, nullptr, 'C'}, + {"reference-mtz", required_argument, nullptr, 'z'}, + {"space-group", required_argument, nullptr, 'S'}, + {"partiality", required_argument, nullptr, 'P'}, + {"anomalous", no_argument, nullptr, 'A'}, + {"refine-bfactor", no_argument, nullptr, 'B'}, + {"wedge", optional_argument, nullptr, 'w'}, + + {"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA}, + {"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD}, + {"spot-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION}, + {"max-spots", required_argument, nullptr, OPT_MAX_SPOTS}, + {"no-beam-center-refine", no_argument, nullptr, OPT_NO_BEAM_CENTER_REFINE}, + {"scale-merge", no_argument, nullptr, OPT_SCALE_MERGE}, + {"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY}, + {"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC}, + {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS}, + + {nullptr, 0, nullptr, 0} +}; + + void trim_in_place(std::string &t) { size_t b = 0; while (b < t.size() && std::isspace(static_cast(t[b]))) b++; @@ -140,6 +185,8 @@ int main(int argc, char **argv) { int nthreads = 1; int start_image = 0; int end_image = -1; // -1 indicates process until end + int image_stride = 1; + bool verbose = false; bool rotation_indexing = false; std::optional rotation_indexing_range; @@ -172,140 +219,132 @@ int main(int argc, char **argv) { } int opt; - while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:MP:AD:C:T:t:Bw::z:p:q:")) != -1) { - switch (opt) { - case 'o': - output_prefix = optarg; - break; - case 'z': - ref_mtz = optarg; - break; - case 'N': - nthreads = atoi(optarg); - break; - case 's': - start_image = atoi(optarg); - break; - case 'e': - end_image = atoi(optarg); - break; - case 'p': - min_partiality = std::stod(optarg); - break; - case 'q': - min_image_cc = std::stod(optarg); - break; - case 'B': - refine_bfactor = true; - break; - case 'w': - refine_wedge = true; - if (optarg) - wedge_for_scaling = std::stod(optarg); - break; - case 'v': - verbose = true; - break; - case 'd': - d_min_spot_finding = atof(optarg); - break; - case 'c': - max_spot_count_override = atoll(optarg); - break; - case 'R': - rotation_indexing = true; - if (optarg) rotation_indexing_range = atof(optarg); - break; - case 'F': - indexing_algorithm = IndexingAlgorithmEnum::FFT; - break; - case 'X': { - std::string alg = optarg ? optarg : ""; - std::transform(alg.begin(), alg.end(), alg.begin(), - [](unsigned char c) { return static_cast(std::tolower(c)); }); + int option_index = 0; + const char *short_opts = "vo:N:s:e:t:R::X:C:z:FABw::S:"; - if (alg == "ffbidx") - indexing_algorithm = IndexingAlgorithmEnum::FFBIDX; - else if (alg == "fft") - indexing_algorithm = IndexingAlgorithmEnum::FFT; - else if (alg == "fftw") - indexing_algorithm = IndexingAlgorithmEnum::FFTW; - else if (alg == "auto") - indexing_algorithm = IndexingAlgorithmEnum::Auto; - else if (alg == "none") - indexing_algorithm = IndexingAlgorithmEnum::None; - else { - logger.Error("Invalid indexing algorithm: {}", alg); - print_usage(logger); - exit(EXIT_FAILURE); - } - break; - } - case 'x': - refine_beam_center = false; - break; - case 'D': - d_min_scale_merge = atof(optarg); - logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_scale_merge.value()); - break; - case 'S': - space_group_number = atoi(optarg); - break; - case 'M': - run_scaling = true; - break; - case 'A': - anomalous_mode = true; - break; - case 'T': - sigma_spot_finding = atof(optarg); - logger.Info("Noise threshold level for spot finding set to {:.2f} sigma", sigma_spot_finding); - break; - case 't': - photon_count_threshold_spot_finding = atoi(optarg); - logger.Info("Photon-count threshold level for spot finding set to {:d}", - photon_count_threshold_spot_finding); - break; - case 'C': { - auto uc = parse_unit_cell_arg(optarg); - if (!uc.has_value()) { - logger.Error( - "Invalid -C unit cell. Expected: -C\"a,b,c,alpha,beta,gamma\" (6 floats, comma-separated, no spaces). Got: {}", - optarg ? optarg : ""); - print_usage(logger); - exit(EXIT_FAILURE); - } - fixed_reference_unit_cell = uc; - logger.Info( - "Fixed reference unit cell set: a={:.3f} b={:.3f} c={:.3f} alpha={:.3f} beta={:.3f} gamma={:.3f}", - uc->a, uc->b, uc->c, uc->alpha, uc->beta, uc->gamma); - break; - } - case 'P': - if (strcmp(optarg, "unity") == 0) - partiality_model = PartialityModel::Unity; - else if (strcmp(optarg, "fixed") == 0) - partiality_model = PartialityModel::Fixed; - else if (strcmp(optarg, "rot") == 0) - partiality_model = PartialityModel::Rotation; - else { - logger.Error("Invalid partiality mode: {}", optarg); - print_usage(logger); - exit(EXIT_FAILURE); - } - break; - case 'i': - scaling_iter = atoi(optarg); - if (scaling_iter <= 0) { - logger.Error("Invalid scaling iteration count: {}", scaling_iter); - exit(EXIT_FAILURE); - } - break; - default: + while ((opt = getopt_long(argc, argv, short_opts, long_options, &option_index)) != -1) { + switch (opt) { + case 'o': + output_prefix = optarg; + break; + case 'v': + verbose = true; + break; + case 'N': + nthreads = atoi(optarg); + break; + case 's': + start_image = atoi(optarg); + break; + case 'e': + end_image = atoi(optarg); + break; + case 't': + image_stride = atoi(optarg); + break; + case 'R': + rotation_indexing = true; + if (optarg) rotation_indexing_range = atof(optarg); + break; + case 'X': { + std::string alg = optarg ? optarg : ""; + std::transform(alg.begin(), alg.end(), alg.begin(), + [](unsigned char c) { return static_cast(std::tolower(c)); }); + + if (alg == "ffbidx") + indexing_algorithm = IndexingAlgorithmEnum::FFBIDX; + else if (alg == "fft") + indexing_algorithm = IndexingAlgorithmEnum::FFT; + else if (alg == "fftw") + indexing_algorithm = IndexingAlgorithmEnum::FFTW; + else if (alg == "auto") + indexing_algorithm = IndexingAlgorithmEnum::Auto; + else if (alg == "none") + indexing_algorithm = IndexingAlgorithmEnum::None; + else { + logger.Error("Invalid indexing algorithm: {}", alg); print_usage(logger); exit(EXIT_FAILURE); + } + break; } + case 'C': { + auto uc = parse_unit_cell_arg(optarg); + if (!uc.has_value()) { + logger.Error( + "Invalid unit cell. Expected: \"a,b,c,alpha,beta,gamma\" (6 floats, comma-separated, no spaces). Got: {}", + optarg ? optarg : ""); + print_usage(logger); + exit(EXIT_FAILURE); + } + fixed_reference_unit_cell = uc; + logger.Info( + "Fixed reference unit cell set: a={:.3f} b={:.3f} c={:.3f} alpha={:.3f} beta={:.3f} gamma={:.3f}", + uc->a, uc->b, uc->c, uc->alpha, uc->beta, uc->gamma); + break; + } + case 'z': + ref_mtz = optarg; + break; + case 'F': + indexing_algorithm = IndexingAlgorithmEnum::FFT; + break; + case 'A': + anomalous_mode = true; + break; + case 'B': + refine_bfactor = true; + break; + case 'w': + refine_wedge = true; + if (optarg) + wedge_for_scaling = std::stod(optarg); + break; + case 'S': + space_group_number = atoi(optarg); + break; + case OPT_SPOT_SIGMA: + sigma_spot_finding = atof(optarg); + logger.Info("Noise threshold level for spot finding set to {:.2f} sigma", sigma_spot_finding); + break; + case OPT_SPOT_THRESHOLD: + photon_count_threshold_spot_finding = atoi(optarg); + logger.Info("Photon-count threshold level for spot finding set to {:d}", + photon_count_threshold_spot_finding); + break; + case OPT_SPOT_RESOLUTION: + d_min_spot_finding = atof(optarg); + logger.Info("High resolution limit for spot finding set to {:.2f} A", d_min_spot_finding); + break; + case OPT_MAX_SPOTS: + max_spot_count_override = atoll(optarg); + logger.Info("Max spot count overridden to {}", max_spot_count_override.value()); + break; + case OPT_NO_BEAM_CENTER_REFINE: + refine_beam_center = false; + break; + case OPT_SCALE_MERGE: + run_scaling = true; + break; + case OPT_MIN_PARTIALITY: + min_partiality = std::stod(optarg); + break; + case OPT_MIN_IMAGE_CC: + min_image_cc = std::stod(optarg); + break; + case OPT_SCALING_ITERATIONS: + scaling_iter = atoi(optarg); + if (scaling_iter <= 0) { + logger.Error("Invalid scaling iteration count: {}", scaling_iter); + exit(EXIT_FAILURE); + } + break; + + default: + print_usage(logger); + exit(EXIT_FAILURE); } +} if (optind != argc - 1) { logger.Error("Input file not specified"); @@ -354,11 +393,16 @@ int main(int argc, char **argv) { if (end_image < 0 || end_image > total_images_in_file) end_image = total_images_in_file; - int images_to_process = end_image - start_image; + if (image_stride < 0) { + logger.Error("Image stride cannot be negative"); + exit(EXIT_FAILURE); + } + + int images_to_process = (end_image - start_image) / image_stride; if (images_to_process <= 0) { - logger.Warning("No images to process (Start: {}, End: {}, Total: {})", start_image, end_image, - total_images_in_file); + logger.Warning("No images to process (Start: {}, End: {} Stride: {}, Total: {})", start_image, end_image, + image_stride, total_images_in_file); return 0; } @@ -483,7 +527,7 @@ int main(int argc, char **argv) { AzimuthalIntegrationProfile profile(mapping); while (true) { - int current_idx_offset = processed_count.fetch_add(1); + int current_idx_offset = processed_count.fetch_add(image_stride); int image_idx = start_image + current_idx_offset; if (image_idx >= end_image) break;