diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 104967bd..2eaac793 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -44,4 +44,4 @@ ADD_SUBDIRECTORY(geom_refinement) ADD_SUBDIRECTORY(lattice_search) ADD_SUBDIRECTORY(scale_merge) -TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement gemmi) +TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge gemmi) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 8c9457ec..fa45a846 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -240,3 +240,32 @@ std::optional IndexAndRefine::Finalize() { return rotation_indexer->GetLattice(); return {}; } + +std::optional IndexAndRefine::ScaleRotationData(const ScaleMergeOptions &opts) const { + std::vector snapshot; + { + std::unique_lock ul(reflections_mutex); + snapshot = reflections; // cheap copy under lock + } + + // Need a reasonable number of reflections to make refinement meaningful + constexpr size_t kMinReflections = 20; + if (snapshot.size() < kMinReflections) + return std::nullopt; + + // Build options focused on mosaicity refinement but allow caller override + ScaleMergeOptions options = opts; + + // If the experiment provides a wedge, propagate it + if (experiment.GetGoniometer().has_value()) + options.wedge_deg = experiment.GetGoniometer()->GetWedge_deg(); + + // If caller left space_group unset, try to pick it from the indexed lattice + if (!options.space_group.has_value()) { + auto sg = experiment.GetGemmiSpaceGroup(); + if (sg) + options.space_group = *sg; + } + + return ScaleAndMergeReflectionsCeres(snapshot, options); +} \ No newline at end of file diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index c3ff50f8..15fb2b13 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -12,6 +12,7 @@ #include "bragg_prediction/BraggPrediction.h" #include "indexing/IndexerThreadPool.h" #include "lattice_search/LatticeSearch.h" +#include "scale_merge/ScaleAndMerge.h" #include "RotationIndexer.h" #include "RotationParameters.h" @@ -43,7 +44,7 @@ class IndexAndRefine { : experiment(experiment_ref) {} }; - std::mutex reflections_mutex; + mutable std::mutex reflections_mutex; std::vector reflections; IndexingOutcome DetermineLatticeAndSymmetry(DataMessage &msg); @@ -56,6 +57,12 @@ class IndexAndRefine { public: IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer); void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction); + + /// Run scale-and-merge on accumulated reflections to refine per-image + /// mosaicity (and optionally B-factors / scale factors). + /// Returns std::nullopt if there are too few reflections to be meaningful. + std::optional ScaleRotationData(const ScaleMergeOptions &opts = {}) const; + std::optional Finalize(); }; diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index acc3f91c..bd60ef07 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include "../reader/JFJochHDF5Reader.h" #include "../common/Logger.h" @@ -27,15 +28,16 @@ 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(" -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(" -R[num] Rotation indexing (optional: min angular range deg)"); + logger.Info(" -R[num] Rotation indexing (optional: min angular range deg)"); logger.Info(" -F Use FFT indexing algorithm (default: Auto)"); logger.Info(" -x No least-square beam center refinement"); - logger.Info(" -d High resolution limit for spot finding (default: 1.5)"); + logger.Info(" -d High resolution limit for spot finding (default: 1.5)"); + logger.Info(" -S Run scaling (refine mosaicity) and write scaled.hkl + image.dat"); } int main(int argc, char **argv) { @@ -55,6 +57,7 @@ int main(int argc, char **argv) { std::optional rotation_indexing_range; bool use_fft = false; bool refine_beam_center = true; + bool run_scaling = false; float d_high = 1.5; @@ -94,6 +97,9 @@ int main(int argc, char **argv) { case 'd': d_high = atof(optarg); break; + case 'S': + run_scaling = true; + break; default: print_usage(logger); exit(EXIT_FAILURE); @@ -354,7 +360,7 @@ int main(int argc, char **argv) { IndexAndRefine global_indexer(experiment, &indexer_pool); const auto rotation_indexer_ret = global_indexer.Finalize(); - if (rotation_indexer_ret.has_value()) { + if (rotation_indexer_ret.has_value()) { end_msg.rotation_lattice = rotation_indexer_ret->lattice; end_msg.rotation_lattice_type = LatticeMessage{ .centering = rotation_indexer_ret->search_result.centering, @@ -364,6 +370,66 @@ int main(int argc, char **argv) { logger.Info("Rotation Indexing found lattice"); } + // --- Optional: run scaling (mosaicity refinement) on accumulated reflections --- + if (run_scaling) { + logger.Info("Running scaling (mosaicity refinement) ..."); + + ScaleMergeOptions scale_opts; + scale_opts.refine_b_factor = false; // B-factor refinement doesn't make sense for rotation + scale_opts.refine_mosaicity = true; + scale_opts.max_num_iterations = 500; + scale_opts.max_solver_time_s = 240.0; // generous cutoff for now + + auto scale_start = std::chrono::steady_clock::now(); + auto scale_result = indexer.ScaleRotationData(scale_opts); + auto scale_end = std::chrono::steady_clock::now(); + double scale_time = std::chrono::duration(scale_end - scale_start).count(); + + if (scale_result) { + logger.Info("Scaling completed in {:.2f} s (GoF² = {:.4f}, {} unique reflections, {} images)", + scale_time, scale_result->gof2, + scale_result->merged.size(), scale_result->image_ids.size()); + + // Write scaled.hkl (h k l I sigma) + { + const std::string hkl_path = output_prefix + "_scaled.hkl"; + std::ofstream hkl_file(hkl_path); + if (!hkl_file) { + logger.Error("Cannot open {} for writing", hkl_path); + } else { + hkl_file << "# h k l I sigma\n"; + for (const auto& r : scale_result->merged) { + hkl_file << r.h << " " << r.k << " " << r.l << " " + << r.I << " " << r.sigma << "\n"; + } + hkl_file.close(); + logger.Info("Wrote {} reflections to {}", scale_result->merged.size(), hkl_path); + } + } + + // Write image.dat (image_id mosaicity_deg K) + { + const std::string img_path = output_prefix + "_image.dat"; + std::ofstream img_file(img_path); + if (!img_file) { + logger.Error("Cannot open {} for writing", img_path); + } else { + img_file << "# image_id mosaicity_deg K\n"; + for (size_t i = 0; i < scale_result->image_ids.size(); ++i) { + img_file << scale_result->image_ids[i] << " " + << scale_result->mosaicity_deg[i] << " " + << scale_result->image_scale_k[i] << "\n"; + } + img_file.close(); + logger.Info("Wrote {} image records to {}", scale_result->image_ids.size(), img_path); + } + } + } else { + logger.Warning("Scaling skipped — too few reflections accumulated (need >= 20)"); + logger.Info("Scaling wall-clock time: {:.2f} s", scale_time); + } + } + // Write End Message writer->WriteHDF5(end_msg); auto stats = writer->Finalize();