diff --git a/broker/OpenAPIConvert.cpp b/broker/OpenAPIConvert.cpp index be884867..61b02aba 100644 --- a/broker/OpenAPIConvert.cpp +++ b/broker/OpenAPIConvert.cpp @@ -1023,6 +1023,10 @@ org::openapitools::server::model::Indexing_settings Convert(const IndexingSettin case GeomRefinementAlgorithmEnum::BeamCenter: refinement.setValue(org::openapitools::server::model::Geom_refinement_algorithm::eGeom_refinement_algorithm::BEAMCENTER); break; + case GeomRefinementAlgorithmEnum::Multi: + // No dedicated API value (offline rugnux only); report it as beam+lattice. + refinement.setValue(org::openapitools::server::model::Geom_refinement_algorithm::eGeom_refinement_algorithm::BEAMCENTER); + break; } ret.setGeomRefinementAlgorithm(refinement); diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index 6e900440..7b92297d 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -6,7 +6,7 @@ #include enum class IndexingAlgorithmEnum {FFBIDX, FFT, FFTW, Auto, None}; -enum class GeomRefinementAlgorithmEnum {None, OrientationOnly, BeamCenter}; +enum class GeomRefinementAlgorithmEnum {None, OrientationOnly, BeamCenter, Multi}; class IndexingSettings { IndexingAlgorithmEnum algorithm; diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 8d542eac..b73a97ab 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -176,6 +176,23 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data return outcome; } +namespace { +// Count spots whose fractional Miller index falls within the indexing tolerance of an integer for a +// given lattice + geometry - the "how well does this model explain the spots" score used by -r multi. +int CountIndexedSpots(const DiffractionGeometry &geom, const CrystalLattice &latt, + const std::vector &spots, float tol_sq) { + const Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2(); + int n = 0; + for (const auto &s : spots) { + const Coord recip = s.ReciprocalCoord(geom); + const float hf = recip * a, kf = recip * b, lf = recip * c; + const float dh = hf - std::round(hf), dk = kf - std::round(kf), dl = lf - std::round(lf); + if (dh * dh + dk * dk + dl * dl < tol_sq) ++n; + } + return n; +} +} // namespace + void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::IndexingOutcome &outcome) { if (!outcome.lattice_candidate) return; @@ -213,6 +230,39 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In outcome.beam_center_updated = true; } break; + case GeomRefinementAlgorithmEnum::Multi: { + // Try all three refinements per image and keep whichever indexes the most spots. Beam+cell + // refinement helps some stills but diverges on sparse spot lists (few spots, long axes), + // where it pushes a good lattice out of tolerance; scoring by indexed-spot count lets each + // image fall back to orientation-only or no refinement when refinement would hurt. Ties + // prefer less refinement (strict >, not >=) to avoid overfitting. + const float tol = experiment.GetIndexingSettings().GetTolerance(); + const float tol_sq = tol * tol; + + XtalOptimizerData d_none = data; + XtalOptimizerData d_orient = data; + XtalOptimizerRotationOnly(d_orient, msg.spots, 0.2); + XtalOptimizerRotationOnly(d_orient, msg.spots, 0.1); + XtalOptimizerRotationOnly(d_orient, msg.spots, 0.05); + XtalOptimizerData d_beam = data; + const bool beam_ok = XtalOptimizer(d_beam, {msg.spots}); + + const int s_none = CountIndexedSpots(d_none.geom, d_none.latt, msg.spots, tol_sq); + const int s_orient = CountIndexedSpots(d_orient.geom, d_orient.latt, msg.spots, tol_sq); + const int s_beam = beam_ok ? CountIndexedSpots(d_beam.geom, d_beam.latt, msg.spots, tol_sq) : -1; + + if (s_beam > s_none && s_beam > s_orient) { + data = d_beam; + outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl()) + .BeamY_pxl(data.geom.GetBeamY_pxl()); + outcome.beam_center_updated = true; + } else if (s_orient > s_none) { + data = d_orient; + } else { + data = d_none; + } + break; + } } outcome.lattice_candidate = data.latt; diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 3ea8f241..b8390bde 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -36,6 +36,7 @@ namespace { switch (r) { case GeomRefinementAlgorithmEnum::None: return "none"; case GeomRefinementAlgorithmEnum::OrientationOnly: return "orientation"; + case GeomRefinementAlgorithmEnum::Multi: return "multi"; case GeomRefinementAlgorithmEnum::BeamCenter: default: return "beam_and_lattice"; } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 42446fec..fa417766 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -75,7 +75,7 @@ void print_usage() { std::cout << " -X, --indexing-algorithm Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << 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 << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice)" << std::endl; + std::cout << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice|multi); multi tries all three per image and keeps whichever indexes the most spots" << std::endl; std::cout << std::endl; std::cout << " Scaling and merging (on by default)" << std::endl; @@ -604,6 +604,8 @@ int main(int argc, char **argv) { refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter; else if (alg == "orientation") refinement_algorithm = GeomRefinementAlgorithmEnum::OrientationOnly; + else if (alg == "multi") + refinement_algorithm = GeomRefinementAlgorithmEnum::Multi; else { logger.Error("Invalid geom refinement algorithm: {}", alg); print_usage(); diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 0497e96e..44ecfeea 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -455,6 +455,9 @@ void NXmx::MX(const StartMessage &start) { case GeomRefinementAlgorithmEnum::BeamCenter: hdf5_file->SaveScalar("/entry/MX/geom_refinement_algorithm", "beam_center"); break; + case GeomRefinementAlgorithmEnum::Multi: + hdf5_file->SaveScalar("/entry/MX/geom_refinement_algorithm", "multi"); + break; default: break; }