v1.0.0-rc.148 (#58)

This is an UNSTABLE release. The release has significant modifications for data processing - in case of troubles go back to 1.0.0-rc.144.

* jfjoch_broker: Improve azimuthal integration (add <I^2> calculation)
* jfjoch_broker: Fixes around indexing, aiming to handle multi-lattice crystals (work in progress, it is not fully integrated)
* jfjoch_writer: Save mean(I), stddev(I), and count(I) for each azimuthal bin

Reviewed-on: #58

tools/jfjoch_process.cpp

@@ -54,6 +54,7 @@ void print_usage() {
     std::cout << "   -R, --two-pass-rotation[=num]    Two-pass offline rotation indexing (optional: number of images, default: 30)" << std::endl;
     std::cout << "   --single-pass-rotation[=num]     Use online-like single-pass rotation indexing (optional: min angular range deg)" << std::endl;
     std::cout << "   --redo-rotation-spots            Redo spot finding for two-pass rotation indexing" << std::endl;
+    std::cout << "   --force-rotation-lattice <vec>   Force rotation indexer with external lattice (in Angstrom) : \"a0x,a0y,a0z,a1x,a1y,a1z,a2x,a2y,a2z\" (9 floats, skips first pass)" << std::endl;
     std::cout << "   -X, --indexing-algorithm <txt>   Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << 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;
@@ -85,7 +86,8 @@ enum {
     OPT_SCALING_HIGH_RESOLUTION,
     OPT_SCALING_OUTPUT,
     OPT_SINGLE_PASS_ROTATION,
-    OPT_REDO_ROTATION_SPOTS
+    OPT_REDO_ROTATION_SPOTS,
+    OPT_FORCE_ROTATION_LATTICE
 };
 
 static option long_options[] = {
@@ -109,6 +111,8 @@ static option long_options[] = {
     {"two-pass-rotation", optional_argument, nullptr, 'R'},
     {"single-pass-rotation", optional_argument, nullptr, OPT_SINGLE_PASS_ROTATION},
     {"redo-rotation-spots", no_argument, nullptr, OPT_REDO_ROTATION_SPOTS},
+    {"force-rotation-lattice", required_argument, nullptr, OPT_FORCE_ROTATION_LATTICE},
+
 
     {"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA},
     {"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD},
@@ -130,6 +134,16 @@ void trim_in_place(std::string &t) {
     t = t.substr(b, e - b);
 };
 
+bool parse_float_strict(const std::string &t, float &out) {
+    try {
+        size_t idx = 0;
+        out = std::stof(t, &idx);
+        return idx == t.size();
+    } catch (...) {
+        return false;
+    }
+};
+
 std::optional<UnitCell> parse_unit_cell_arg(const char *arg) {
     if (!arg)
         return std::nullopt;
@@ -160,15 +174,7 @@ std::optional<UnitCell> parse_unit_cell_arg(const char *arg) {
     if (parts.size() != 6)
         return std::nullopt;
 
-    auto parse_float_strict = [](const std::string &t, float &out) -> bool {
-        try {
-            size_t idx = 0;
-            out = std::stof(t, &idx);
-            return idx == t.size();
-        } catch (...) {
-            return false;
-        }
-    };
+
 
     UnitCell uc{};
     if (!parse_float_strict(parts[0], uc.a)) return std::nullopt;
@@ -181,6 +187,43 @@ std::optional<UnitCell> parse_unit_cell_arg(const char *arg) {
     return uc;
 }
 
+std::optional<CrystalLattice> parse_lattice_arg(const char *arg) {
+    if (!arg)
+        return std::nullopt;
+
+    std::string s(arg);
+    trim_in_place(s);
+
+    if (s.size() >= 2 && ((s.front() == '"' && s.back() == '"') || (s.front() == '\'' && s.back() == '\''))) {
+        s = s.substr(1, s.size() - 2);
+        trim_in_place(s);
+    }
+
+    std::vector<std::string> parts;
+    parts.reserve(9);
+    size_t start = 0;
+    while (true) {
+        size_t pos = s.find(',', start);
+        if (pos == std::string::npos) {
+            parts.push_back(s.substr(start));
+            break;
+        }
+        parts.push_back(s.substr(start, pos - start));
+        start = pos + 1;
+    }
+
+    if (parts.size() != 9)
+        return std::nullopt;
+
+    std::vector<float> vals(9);
+    for (int i = 0; i < 9; i++) {
+        if (!parse_float_strict(parts[i], vals[i]))
+            return std::nullopt;
+    }
+
+    return CrystalLattice(vals);
+}
+
 std::vector<int> select_equally_spaced_image_ordinals(int images_to_process, int requested_images) {
     std::vector<int> ret;
 
@@ -248,6 +291,7 @@ int main(int argc, char **argv) {
     double min_partiality = 0.02;
     double min_image_cc = 0.0;
     int64_t scaling_iter = 3;
+    std::optional<CrystalLattice> forced_rotation_lattice;
 
     IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto;
     GeomRefinementAlgorithmEnum refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter;
@@ -312,6 +356,28 @@ int main(int argc, char **argv) {
         case OPT_REDO_ROTATION_SPOTS:
             reuse_rotation_spots = false;
             break;
+        case OPT_FORCE_ROTATION_LATTICE: {
+            if (rotation_indexing) {
+                logger.Error("Rotation indexing already enabled");
+                exit(EXIT_FAILURE);
+            }
+            rotation_indexing = true;
+
+            auto latt = parse_lattice_arg(optarg);
+            if (!latt.has_value()) {
+                logger.Error(
+                    "Invalid rotation lattice. Expected: \"a0x,a0y,a0z,a1x,a1y,a1z,a2x,a2y,a2z\" (9 floats, comma-separated). Got: {}",
+                    optarg ? optarg : "<null>");
+                print_usage();
+                exit(EXIT_FAILURE);
+            }
+            forced_rotation_lattice = latt;
+            auto uc = latt->GetUnitCell();
+            logger.Info(
+                "Forced rotation lattice 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 'X': {
             std::string alg = optarg ? optarg : "";
             std::transform(alg.begin(), alg.end(), alg.begin(),
@@ -635,7 +701,10 @@ int main(int argc, char **argv) {
     if (!reference_data.empty())
         indexer.ReferenceIntensities(reference_data);
 
-    if (rotation_indexing && two_pass_rotation) {
+    if (forced_rotation_lattice.has_value()) {
+        indexer.ForceRotationIndexerLattice(*forced_rotation_lattice);
+        logger.Info("Rotation indexer lattice forced externally - skipping first pass indexing");
+    } else if (rotation_indexing && two_pass_rotation) {
         const auto selected_ordinals = select_equally_spaced_image_ordinals(
             images_to_process, rotation_indexing_image_count);