275 lines
10 KiB
C++
275 lines
10 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "IndexAndRefine.h"
|
|
|
|
#include "bragg_integration/BraggIntegrate2D.h"
|
|
#include "bragg_integration/CalcISigma.h"
|
|
#include "geom_refinement/XtalOptimizer.h"
|
|
#include "indexing/AnalyzeIndexing.h"
|
|
#include "indexing/FFTIndexer.h"
|
|
#include "lattice_search/LatticeSearch.h"
|
|
|
|
IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer)
|
|
: index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()),
|
|
experiment(x),
|
|
geom_(x.GetDiffractionGeometry()),
|
|
indexer_(indexer) {
|
|
if (indexer && x.IsRotationIndexing())
|
|
rotation_indexer = std::make_unique<RotationIndexer>(x, *indexer);
|
|
reflections.resize(x.GetImageNum());
|
|
mosaicity.resize(x.GetImageNum(), NAN);
|
|
}
|
|
|
|
IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(DataMessage &msg) {
|
|
IndexingOutcome outcome(experiment);
|
|
|
|
if (rotation_indexer) {
|
|
auto result = rotation_indexer->ProcessImage(msg.number, msg.spots);
|
|
if (result) {
|
|
|
|
// get rotated lattice
|
|
auto gon = result->axis;
|
|
if (gon) {
|
|
const float angle_deg = gon->GetAngle_deg(msg.number) + gon->GetWedge_deg() / 2.0f;
|
|
outcome.lattice_candidate = result->lattice.Multiply(gon->GetTransformationAngle(-angle_deg));
|
|
}
|
|
|
|
outcome.experiment.BeamX_pxl(result->geom.GetBeamX_pxl())
|
|
.BeamY_pxl(result->geom.GetBeamY_pxl())
|
|
.DetectorDistance_mm(result->geom.GetDetectorDistance_mm())
|
|
.PoniRot1_rad(result->geom.GetPoniRot1_rad())
|
|
.PoniRot2_rad(result->geom.GetPoniRot2_rad())
|
|
.Goniometer(result->axis);
|
|
outcome.symmetry.centering = result->search_result.centering;
|
|
outcome.symmetry.niggli_class = result->search_result.niggli_class;
|
|
outcome.symmetry.crystal_system = result->search_result.system;
|
|
}
|
|
return outcome;
|
|
}
|
|
|
|
// Convert input spots to reciprocal space
|
|
std::vector<Coord> recip;
|
|
recip.reserve(msg.spots.size());
|
|
for (const auto &i: msg.spots) {
|
|
if (index_ice_rings || !i.ice_ring)
|
|
recip.push_back(i.ReciprocalCoord(geom_));
|
|
}
|
|
|
|
auto indexer_result = indexer_->Run(experiment, recip).get();
|
|
msg.indexing_time_s = indexer_result.indexing_time_s;
|
|
|
|
if (indexer_result.lattice.empty())
|
|
return outcome;
|
|
|
|
auto latt = indexer_result.lattice[0];
|
|
auto sg = experiment.GetGemmiSpaceGroup();
|
|
|
|
// If space group provided => always enforce symmetry in refinement
|
|
// If space group not provided => guess symmetry
|
|
if (sg) {
|
|
// If space group provided but no unit cell fixed, it is better to keep triclinic for now
|
|
if (experiment.GetUnitCell()) {
|
|
outcome.symmetry = LatticeMessage{
|
|
.centering = sg->centring_type(),
|
|
.niggli_class = 0,
|
|
.crystal_system = sg->crystal_system()
|
|
};
|
|
}
|
|
outcome.lattice_candidate = latt;
|
|
} else {
|
|
auto sym_result = LatticeSearch(latt);
|
|
outcome.symmetry = LatticeMessage{
|
|
.centering = sym_result.centering,
|
|
.niggli_class = sym_result.niggli_class,
|
|
.crystal_system = sym_result.system
|
|
};
|
|
outcome.lattice_candidate = sym_result.conventional;
|
|
}
|
|
|
|
return outcome;
|
|
}
|
|
|
|
void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::IndexingOutcome &outcome) {
|
|
if (!outcome.lattice_candidate)
|
|
return;
|
|
|
|
XtalOptimizerData data{
|
|
.geom = outcome.experiment.GetDiffractionGeometry(),
|
|
.latt = *outcome.lattice_candidate,
|
|
.crystal_system = outcome.symmetry.crystal_system,
|
|
.min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(),
|
|
.refine_beam_center = true,
|
|
.refine_distance_mm = false,
|
|
.refine_detector_angles = false,
|
|
.max_time = 0.04 // 40 ms is max allowed time for the operation
|
|
};
|
|
|
|
if (experiment.IsRotationIndexing()) {
|
|
data.refine_beam_center = false;
|
|
data.refine_rotation_axis = false;
|
|
data.refine_unit_cell = false;
|
|
}
|
|
|
|
if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Trigonal)
|
|
data.crystal_system = gemmi::CrystalSystem::Hexagonal;
|
|
|
|
switch (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm()) {
|
|
case GeomRefinementAlgorithmEnum::None:
|
|
break;
|
|
case GeomRefinementAlgorithmEnum::BeamCenter:
|
|
if (XtalOptimizer(data, msg.spots)) {
|
|
outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl())
|
|
.BeamY_pxl(data.geom.GetBeamY_pxl());
|
|
outcome.beam_center_updated = true;
|
|
}
|
|
break;
|
|
}
|
|
|
|
outcome.lattice_candidate = data.latt;
|
|
|
|
if (outcome.beam_center_updated) {
|
|
msg.beam_corr_x = data.beam_corr_x;
|
|
msg.beam_corr_y = data.beam_corr_y;
|
|
}
|
|
}
|
|
|
|
void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
|
|
const SpotFindingSettings &spot_finding_settings,
|
|
const CompressedImage &image,
|
|
BraggPrediction &prediction,
|
|
const IndexAndRefine::IndexingOutcome &outcome) {
|
|
if (!spot_finding_settings.quick_integration)
|
|
return;
|
|
if (!outcome.lattice_candidate)
|
|
return;
|
|
|
|
CrystalLattice latt = outcome.lattice_candidate.value();
|
|
|
|
if (rotation_indexer) {
|
|
// Use moving average for mosaicity and profile_radius (also add beam center later)
|
|
if (msg.mosaicity_deg)
|
|
msg.mosaicity_deg = rotation_parameters.Mosaicity(msg.mosaicity_deg.value());
|
|
if (msg.profile_radius) {
|
|
msg.profile_radius = rotation_parameters.ProfileRadius(msg.profile_radius.value());
|
|
}
|
|
}
|
|
|
|
float ewald_dist_cutoff = 0.001f;
|
|
if (msg.profile_radius)
|
|
ewald_dist_cutoff = msg.profile_radius.value() * 2.0f;
|
|
|
|
if (experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA())
|
|
ewald_dist_cutoff = experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA().value() * 3.0f;
|
|
|
|
float wedge_deg = 0.0f;
|
|
float mos_deg = 0.1f;
|
|
|
|
if (experiment.GetGoniometer().has_value()) {
|
|
wedge_deg = experiment.GetGoniometer()->GetWedge_deg() / 2.0;
|
|
|
|
if (msg.mosaicity_deg) {
|
|
mos_deg = msg.mosaicity_deg.value();
|
|
mosaicity[msg.number] = mos_deg;
|
|
}
|
|
}
|
|
|
|
const BraggPredictionSettings settings_prediction{
|
|
.high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
|
|
.ewald_dist_cutoff = ewald_dist_cutoff,
|
|
.max_hkl = 100,
|
|
.centering = outcome.symmetry.centering,
|
|
.wedge_deg = std::fabs(wedge_deg),
|
|
.mosaicity_deg = std::fabs(mos_deg)
|
|
};
|
|
|
|
auto nrefl = prediction.Calc(outcome.experiment, latt, settings_prediction);
|
|
auto refl_ret = BraggIntegrate2D(outcome.experiment, image, prediction.GetReflections(), nrefl, msg.number);
|
|
|
|
constexpr size_t kMaxReflections = 10000;
|
|
if (refl_ret.size() > kMaxReflections) {
|
|
// Keep only smallest d (highest resolution)
|
|
std::nth_element(refl_ret.begin(),
|
|
refl_ret.begin() + static_cast<long>(kMaxReflections),
|
|
refl_ret.end(),
|
|
[](const Reflection& a, const Reflection& b) {
|
|
return a.d < b.d;
|
|
});
|
|
refl_ret.resize(kMaxReflections);
|
|
|
|
// Optional: make output ordered by d (nice for downstream / debugging)
|
|
std::sort(refl_ret.begin(), refl_ret.end(),
|
|
[](const Reflection& a, const Reflection& b) { return a.d < b.d; });
|
|
}
|
|
|
|
{
|
|
std::unique_lock ul(reflections_mutex);
|
|
reflections[msg.number] = refl_ret; // Image is not processed twice, so thread-safe in principle, but better safe than sorry :)
|
|
}
|
|
|
|
msg.reflections = std::move(refl_ret);
|
|
|
|
CalcISigma(msg);
|
|
CalcWilsonBFactor(msg);
|
|
}
|
|
|
|
void IndexAndRefine::ProcessImage(DataMessage &msg,
|
|
const SpotFindingSettings &spot_finding_settings,
|
|
const CompressedImage &image,
|
|
BraggPrediction &prediction) {
|
|
if (!indexer_ || !spot_finding_settings.indexing)
|
|
return;
|
|
|
|
msg.indexing_result = false;
|
|
|
|
IndexingOutcome outcome = DetermineLatticeAndSymmetry(msg);
|
|
if (!outcome.lattice_candidate)
|
|
return;
|
|
|
|
RefineGeometryIfNeeded(msg, outcome);
|
|
if (!outcome.lattice_candidate)
|
|
return;
|
|
|
|
if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate))
|
|
return;
|
|
|
|
msg.lattice_type = outcome.symmetry;
|
|
|
|
QuickPredictAndIntegrate(msg, spot_finding_settings, image, prediction, outcome);
|
|
}
|
|
|
|
std::optional<RotationIndexerResult> IndexAndRefine::Finalize() {
|
|
if (rotation_indexer)
|
|
return rotation_indexer->GetLattice();
|
|
return {};
|
|
}
|
|
|
|
std::optional<ScaleMergeResult> IndexAndRefine::ScaleRotationData(const ScaleMergeOptions &opts) const {
|
|
size_t nrefl = 0;
|
|
for (const auto &i: reflections)
|
|
nrefl += i.size();
|
|
|
|
// Need a reasonable number of reflections to make refinement meaningful
|
|
constexpr size_t kMinReflections = 20;
|
|
if (nrefl < 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();
|
|
options.mosaicity_init_deg_vec = mosaicity;
|
|
}
|
|
|
|
// 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(reflections, options);
|
|
}
|