jfjoch_process: Introduce long options + image stride
This commit is contained in:
+210
-166
@@ -11,6 +11,7 @@
|
||||
#include <chrono>
|
||||
#include <fstream>
|
||||
#include <sstream>
|
||||
#include <getopt.h>
|
||||
|
||||
#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 {<options>} <input.h5>");
|
||||
logger.Info("Options:");
|
||||
logger.Info(" -o<txt> Output file prefix (default: output)");
|
||||
logger.Info(" -N<num> Number of threads (default: 1)");
|
||||
logger.Info(" -s<num> Start image number (default: 0)");
|
||||
logger.Info(" -e<num> End image number (default: all)");
|
||||
logger.Info(" -v Verbose output");
|
||||
logger.Info("");
|
||||
std::cout << "Usage ./jfjoch_analysis {<options>} <input.h5>" << std::endl;
|
||||
std::cout << "Options:" << std::endl;
|
||||
std::cout << " -o, --output-prefix <txt> Output file prefix (default: output)" << std::endl;
|
||||
std::cout << " -N, --threads <num> Number of threads (default: 1)" << std::endl;
|
||||
std::cout << " -s, --start-image <num> Start image number (default: 0)" << std::endl;
|
||||
std::cout << " -e, --end-image <num> End image number (default: all)" << std::endl;
|
||||
std::cout << " -t, --stride <num> 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<num> Noise sigma level for spot finding (default: 3.0)");
|
||||
logger.Info(" -t<num> Photon count threshold for spot finding (default: 10)");
|
||||
logger.Info(" -d<num> High resolution limit for spot finding (default: 1.5)");
|
||||
logger.Info(" -c<num> Max spot count (default: 250)");
|
||||
logger.Info("");
|
||||
logger.Info(" Indexing");
|
||||
logger.Info(" -R[num] Rotation indexing (optional: min angular range deg)");
|
||||
logger.Info(" -X<txt> Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)");
|
||||
logger.Info(" -F Use FFT indexing algorithm (shortcut for -XFFT)");
|
||||
logger.Info(" -S<num> Space group number - used for both indexing and scaling");
|
||||
logger.Info(
|
||||
" -C<cell> 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 <num> Noise sigma level for spot finding (default: 3.0)" << std::endl;
|
||||
std::cout << " --spot-threshold <num> Photon count threshold for spot finding (default: 10)" << std::endl;
|
||||
std::cout << " --spot-resolution <num> High resolution limit for spot finding (default: 1.5)" << std::endl;
|
||||
std::cout << " --max-spots <num> Max spot count (default: 250)" << std::endl;
|
||||
std::cout << std::endl;
|
||||
|
||||
logger.Info("");
|
||||
logger.Info(" Scaling and merging");
|
||||
logger.Info(" -D<num> 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<txt> 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<num> Refine image wedge during scaling with starting wedge value");
|
||||
logger.Info(" -p<num> Minimum partiality to accept reflection (default: 0.02)");
|
||||
logger.Info(" -q<num> Per-image CC limit in percent (default: no limit)");
|
||||
logger.Info(" -i<num> 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 <txt> 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 <num> Space group number - used for both indexing and scaling" << std::endl;
|
||||
std::cout << " -C, --unit-cell <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 <txt> 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 <num> Minimum partiality to accept reflection (default: 0.02)" << std::endl;
|
||||
std::cout << " --min-image-cc <num> Per-image CC limit in percent (default: no limit)" << std::endl;
|
||||
std::cout << " --scaling-iterations <num> Number of scaling iterations with no reference data (default: 3)" << std::endl;
|
||||
std::cout << " -z, --reference-mtz <file> 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<unsigned char>(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<float> 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<char>(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 : "<null>");
|
||||
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<char>(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 : "<null>");
|
||||
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;
|
||||
|
||||
Reference in New Issue
Block a user