ScaleAndMerge: clean-up (no log residual, no fixed mosaicity, return per-image values)
This commit is contained in:
@@ -17,6 +17,7 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool
|
||||
indexer_(indexer) {
|
||||
if (indexer && x.IsRotationIndexing())
|
||||
rotation_indexer = std::make_unique<RotationIndexer>(x, *indexer);
|
||||
reflections.resize(x.GetImageNum());
|
||||
}
|
||||
|
||||
IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(DataMessage &msg) {
|
||||
@@ -198,16 +199,15 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
|
||||
[](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);
|
||||
|
||||
// Append reflections to the class-wide reflections buffer (thread-safe)
|
||||
{
|
||||
std::unique_lock ul(reflections_mutex);
|
||||
reflections.insert(reflections.end(), msg.reflections.begin(), msg.reflections.end());
|
||||
}
|
||||
}
|
||||
|
||||
void IndexAndRefine::ProcessImage(DataMessage &msg,
|
||||
@@ -242,15 +242,13 @@ std::optional<RotationIndexerResult> IndexAndRefine::Finalize() {
|
||||
}
|
||||
|
||||
std::optional<ScaleMergeResult> IndexAndRefine::ScaleRotationData(const ScaleMergeOptions &opts) const {
|
||||
std::vector<Reflection> snapshot;
|
||||
{
|
||||
std::unique_lock ul(reflections_mutex);
|
||||
snapshot = reflections; // cheap copy under lock
|
||||
}
|
||||
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 (snapshot.size() < kMinReflections)
|
||||
if (nrefl < kMinReflections)
|
||||
return std::nullopt;
|
||||
|
||||
// Build options focused on mosaicity refinement but allow caller override
|
||||
@@ -267,5 +265,5 @@ std::optional<ScaleMergeResult> IndexAndRefine::ScaleRotationData(const ScaleMer
|
||||
options.space_group = *sg;
|
||||
}
|
||||
|
||||
return ScaleAndMergeReflectionsCeres(snapshot, options);
|
||||
}
|
||||
return ScaleAndMergeReflectionsCeres(reflections, options);
|
||||
}
|
||||
|
||||
@@ -45,7 +45,7 @@ class IndexAndRefine {
|
||||
};
|
||||
|
||||
mutable std::mutex reflections_mutex;
|
||||
std::vector<Reflection> reflections;
|
||||
std::vector<std::vector<Reflection>> reflections;
|
||||
|
||||
IndexingOutcome DetermineLatticeAndSymmetry(DataMessage &msg);
|
||||
void RefineGeometryIfNeeded(DataMessage &msg, IndexingOutcome &outcome);
|
||||
|
||||
@@ -187,98 +187,86 @@ namespace {
|
||||
};
|
||||
} // namespace
|
||||
|
||||
ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &observations,
|
||||
ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Reflection>> &observations,
|
||||
const ScaleMergeOptions &opt) {
|
||||
ScaleMergeResult out;
|
||||
const bool refine_partiality = opt.wedge_deg > 0.0;
|
||||
ceres::Problem problem;
|
||||
|
||||
struct ObsRef {
|
||||
const Reflection *r = nullptr;
|
||||
int img_id = 0;
|
||||
int img_slot = -1;
|
||||
int hkl_slot = -1;
|
||||
double sigma = 0.0;
|
||||
};
|
||||
|
||||
std::vector<ObsRef> obs;
|
||||
obs.reserve(observations.size());
|
||||
size_t nrefl = 0;
|
||||
for (const auto &i: observations)
|
||||
nrefl += i.size();
|
||||
|
||||
std::unordered_map<int, int> imgIdToSlot;
|
||||
imgIdToSlot.reserve(256);
|
||||
std::vector<ObsRef> obs;
|
||||
obs.reserve(nrefl);
|
||||
|
||||
std::unordered_map<HKLKey, int, HKLKeyHash> hklToSlot;
|
||||
hklToSlot.reserve(observations.size());
|
||||
hklToSlot.reserve(nrefl);
|
||||
|
||||
for (const auto &r: observations) {
|
||||
const double d = SafeD(r.d);
|
||||
if (!std::isfinite(d))
|
||||
continue;
|
||||
if (!std::isfinite(r.I))
|
||||
continue;
|
||||
std::vector<double> g(observations.size(), 1.0);
|
||||
std::vector<double> mosaicity(observations.size(), opt.mosaicity_init_deg);
|
||||
std::vector<uint8_t> image_slot_used(observations.size(), 0);
|
||||
|
||||
if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A)
|
||||
continue;
|
||||
for (int i = 0; i < observations.size(); i++) {
|
||||
for (const auto &r: observations[i]) {
|
||||
const double d = SafeD(r.d);
|
||||
if (!std::isfinite(d))
|
||||
continue;
|
||||
if (!std::isfinite(r.I))
|
||||
continue;
|
||||
|
||||
if (!std::isfinite(r.zeta) || r.zeta <= 0.0f)
|
||||
continue;
|
||||
if (!std::isfinite(r.rlp) || r.rlp == 0.0f)
|
||||
continue;
|
||||
if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A)
|
||||
continue;
|
||||
|
||||
const double sigma = SafeSigma(static_cast<double>(r.sigma), opt.min_sigma);
|
||||
if (!std::isfinite(r.zeta) || r.zeta <= 0.0f)
|
||||
continue;
|
||||
if (!std::isfinite(r.rlp) || r.rlp == 0.0f)
|
||||
continue;
|
||||
|
||||
const int img_id = RoundImageId(r.image_number, opt.image_number_rounding);
|
||||
const double sigma = SafeSigma(r.sigma, opt.min_sigma);
|
||||
|
||||
int img_slot;
|
||||
{
|
||||
auto it = imgIdToSlot.find(img_id);
|
||||
if (it == imgIdToSlot.end()) {
|
||||
img_slot = static_cast<int>(imgIdToSlot.size());
|
||||
imgIdToSlot.emplace(img_id, img_slot);
|
||||
} else {
|
||||
img_slot = it->second;
|
||||
const int img_id = i;
|
||||
|
||||
image_slot_used[img_id] = 1;
|
||||
|
||||
int hkl_slot;
|
||||
try {
|
||||
const HKLKey key = CanonicalizeHKLKey(r, opt);
|
||||
auto it = hklToSlot.find(key);
|
||||
if (it == hklToSlot.end()) {
|
||||
hkl_slot = static_cast<int>(hklToSlot.size());
|
||||
hklToSlot.emplace(key, hkl_slot);
|
||||
} else {
|
||||
hkl_slot = it->second;
|
||||
}
|
||||
} catch (...) {
|
||||
continue;
|
||||
}
|
||||
}
|
||||
|
||||
int hkl_slot;
|
||||
try {
|
||||
const HKLKey key = CanonicalizeHKLKey(r, opt);
|
||||
auto it = hklToSlot.find(key);
|
||||
if (it == hklToSlot.end()) {
|
||||
hkl_slot = static_cast<int>(hklToSlot.size());
|
||||
hklToSlot.emplace(key, hkl_slot);
|
||||
} else {
|
||||
hkl_slot = it->second;
|
||||
}
|
||||
} catch (...) {
|
||||
continue;
|
||||
ObsRef o;
|
||||
o.r = &r;
|
||||
o.img_id = img_id;
|
||||
o.hkl_slot = hkl_slot;
|
||||
o.sigma = sigma;
|
||||
obs.push_back(o);
|
||||
}
|
||||
|
||||
ObsRef o;
|
||||
o.r = &r;
|
||||
o.img_id = img_id;
|
||||
o.img_slot = img_slot;
|
||||
o.hkl_slot = hkl_slot;
|
||||
o.sigma = sigma;
|
||||
obs.push_back(o);
|
||||
}
|
||||
|
||||
const int nimg = static_cast<int>(imgIdToSlot.size());
|
||||
const int nhkl = static_cast<int>(hklToSlot.size());
|
||||
|
||||
out.image_scale_g.assign(nimg, 1.0);
|
||||
out.image_ids.assign(nimg, 0);
|
||||
|
||||
for (const auto &kv: imgIdToSlot) {
|
||||
out.image_ids[kv.second] = kv.first;
|
||||
}
|
||||
|
||||
std::vector<double> g(nimg, 1.0);
|
||||
std::vector<double> Itrue(nhkl, 0.0);
|
||||
|
||||
// Mosaicity: always per-image
|
||||
std::vector<double> mosaicity(nimg, opt.mosaicity_init_deg);
|
||||
for (int i = 0; i < nimg && i < opt.mosaicity_init_deg_vec.size(); ++i) {
|
||||
if (opt.mosaicity_init_deg_vec[i] > 0.0)
|
||||
mosaicity[i] = opt.mosaicity_init_deg_vec[i];
|
||||
for (int i = 0; i < observations.size(); i++) {
|
||||
if (!image_slot_used[i]) {
|
||||
mosaicity[i] = NAN;
|
||||
g[i] = NAN;
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize Itrue from per-HKL median of observed intensities
|
||||
@@ -301,61 +289,73 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &ob
|
||||
}
|
||||
}
|
||||
|
||||
ceres::Problem problem;
|
||||
|
||||
const bool refine_partiality = opt.wedge_deg > 0.0;
|
||||
|
||||
std::vector<bool> is_valid_hkl_slot(nhkl, false);
|
||||
|
||||
for (const auto &o: obs) {
|
||||
size_t mos_slot = opt.per_image_mosaicity ? o.img_slot : 0;
|
||||
auto *cost = new ceres::AutoDiffCostFunction<IntensityResidual, 1, 1, 1, 1>(
|
||||
new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality));
|
||||
|
||||
problem.AddResidualBlock(cost,
|
||||
nullptr,
|
||||
&g[o.img_slot],
|
||||
&mosaicity[mos_slot],
|
||||
&g[o.img_id],
|
||||
&mosaicity[o.img_id],
|
||||
&Itrue[o.hkl_slot]);
|
||||
is_valid_hkl_slot[o.hkl_slot] = true;
|
||||
}
|
||||
|
||||
if (opt.smoothen_g) {
|
||||
for (int i = 0; i < nimg - 2; ++i) {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<SmoothnessRegularizationResidual, 1, 1, 1, 1>(
|
||||
new SmoothnessRegularizationResidual(0.05));
|
||||
|
||||
problem.AddResidualBlock(cost, nullptr,
|
||||
&g[i],
|
||||
&g[i + 1],
|
||||
&g[i + 2]);
|
||||
for (int i = 0; i < g.size(); ++i) {
|
||||
if (image_slot_used[i]) {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<ScaleRegularizationResidual, 1, 1>(
|
||||
new ScaleRegularizationResidual(0.05));
|
||||
problem.AddResidualBlock(cost, nullptr, &g[i]);
|
||||
}
|
||||
}
|
||||
|
||||
if (opt.smoothen_mos && opt.per_image_mosaicity) {
|
||||
for (int i = 0; i < nimg - 2; ++i) {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<SmoothnessRegularizationResidual, 1, 1, 1, 1>(
|
||||
new SmoothnessRegularizationResidual(0.05));
|
||||
|
||||
problem.AddResidualBlock(cost, nullptr,
|
||||
&mosaicity[i],
|
||||
&mosaicity[i + 1],
|
||||
&mosaicity[i + 2]);
|
||||
if (opt.smoothen_g) {
|
||||
for (int i = 0; i < g.size() - 2; ++i) {
|
||||
if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<SmoothnessRegularizationResidual, 1, 1, 1, 1>(
|
||||
new SmoothnessRegularizationResidual(0.05));
|
||||
|
||||
problem.AddResidualBlock(cost, nullptr,
|
||||
&g[i],
|
||||
&g[i + 1],
|
||||
&g[i + 2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (opt.smoothen_mos && opt.refine_mosaicity) {
|
||||
for (int i = 0; i < mosaicity.size() - 2; ++i) {
|
||||
if (image_slot_used[i] && image_slot_used[i + 1] && image_slot_used[i + 2]) {
|
||||
auto *cost = new ceres::AutoDiffCostFunction<SmoothnessRegularizationResidual, 1, 1, 1, 1>(
|
||||
new SmoothnessRegularizationResidual(0.05));
|
||||
|
||||
problem.AddResidualBlock(cost, nullptr,
|
||||
&mosaicity[i],
|
||||
&mosaicity[i + 1],
|
||||
&mosaicity[i + 2]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Scaling factors must be always positive
|
||||
for (int i = 0; i < nimg; ++i)
|
||||
problem.SetParameterLowerBound(&g[i], 0, 1e-12);
|
||||
for (int i = 0; i < g.size(); i++) {
|
||||
if (image_slot_used[i])
|
||||
problem.SetParameterLowerBound(&g[i], 0, 1e-12);
|
||||
}
|
||||
|
||||
// Mosaicity refinement + bounds
|
||||
if (!opt.refine_mosaicity) {
|
||||
for (int i = 0; i < (opt.per_image_mosaicity ? nimg : 1); ++i)
|
||||
problem.SetParameterBlockConstant(&mosaicity[i]);
|
||||
for (int i = 0; i < mosaicity.size(); ++i) {
|
||||
if (image_slot_used[i])
|
||||
problem.SetParameterBlockConstant(&mosaicity[i]);
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < (opt.per_image_mosaicity ? nimg : 1); ++i) {
|
||||
problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg);
|
||||
problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg);
|
||||
for (int i = 0; i < mosaicity.size(); ++i) {
|
||||
if (image_slot_used[i]) {
|
||||
problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg);
|
||||
problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -375,18 +375,15 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &ob
|
||||
ceres::Solver::Summary summary;
|
||||
ceres::Solve(options, &problem, &summary);
|
||||
|
||||
// --- Export per-image results ---
|
||||
for (int i = 0; i < nimg; ++i)
|
||||
out.image_scale_g[i] = g[i];
|
||||
|
||||
out.mosaicity_deg.resize(nimg);
|
||||
for (int i = 0; i < nimg; ++i)
|
||||
out.mosaicity_deg[i] = opt.per_image_mosaicity ? mosaicity[i] : mosaicity[0];
|
||||
ScaleMergeResult out;
|
||||
out.image_scale_g = g;
|
||||
out.mosaicity_deg.resize(observations.size());
|
||||
for (int i = 0; i < out.mosaicity_deg.size(); i++)
|
||||
out.mosaicity_deg[i] = mosaicity[i];
|
||||
|
||||
std::vector<HKLKey> slotToHKL(nhkl);
|
||||
for (const auto &kv: hklToSlot) {
|
||||
for (const auto &kv: hklToSlot)
|
||||
slotToHKL[kv.second] = kv.first;
|
||||
}
|
||||
|
||||
out.merged.resize(nhkl);
|
||||
for (int h = 0; h < nhkl; ++h) {
|
||||
@@ -432,15 +429,14 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &ob
|
||||
for (const auto &o: obs) {
|
||||
const Reflection &r = *o.r;
|
||||
const double lp = SafeInv(static_cast<double>(r.rlp), 1.0);
|
||||
const double G_i = g[o.img_slot];
|
||||
const double G_i = g[o.img_id];
|
||||
|
||||
// Compute partiality with refined mosaicity
|
||||
double partiality;
|
||||
const size_t mos_slot = opt.per_image_mosaicity ? o.img_slot : 0;
|
||||
if (refine_partiality && mosaicity[mos_slot] > 0.0) {
|
||||
if (refine_partiality && mosaicity[o.img_id] > 0.0) {
|
||||
const double c1 = r.zeta / std::sqrt(2.0);
|
||||
const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[mos_slot];
|
||||
const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[mos_slot];
|
||||
const double arg_plus = (r.delta_phi_deg + half_wedge) * c1 / mosaicity[o.img_id];
|
||||
const double arg_minus = (r.delta_phi_deg - half_wedge) * c1 / mosaicity[o.img_id];
|
||||
partiality = (std::erf(arg_plus) - std::erf(arg_minus)) / 2.0;
|
||||
} else {
|
||||
partiality = r.partiality;
|
||||
|
||||
@@ -32,7 +32,6 @@ struct ScaleMergeOptions {
|
||||
|
||||
// --- Mosaicity (user input in degrees; internally converted to radians) ---
|
||||
bool refine_mosaicity = true;
|
||||
bool per_image_mosaicity = true;
|
||||
|
||||
double mosaicity_init_deg = 0.17; // ~0.003 rad
|
||||
double mosaicity_min_deg = 1e-3;
|
||||
@@ -81,14 +80,11 @@ struct MergeStatistics {
|
||||
struct ScaleMergeResult {
|
||||
std::vector<MergedReflection> merged;
|
||||
std::vector<double> image_scale_g;
|
||||
std::vector<int> image_ids;
|
||||
|
||||
// One mosaicity value per image (degrees).
|
||||
std::vector<double> mosaicity_deg;
|
||||
|
||||
/// Per-shell and overall merging statistics (populated after merging)
|
||||
MergeStatistics statistics;
|
||||
};
|
||||
|
||||
ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& observations,
|
||||
ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Reflection>>& observations,
|
||||
const ScaleMergeOptions& opt = {});
|
||||
@@ -42,7 +42,7 @@ void print_usage(Logger &logger) {
|
||||
logger.Info(" -D<num> High resolution limit for scaling/merging (default: 0.0; no limit)");
|
||||
logger.Info(" -S<num> Space group number");
|
||||
logger.Info(" -M Scale and merge (refine mosaicity) and write scaled.hkl + image.dat");
|
||||
logger.Info(" -m<txt> Mosaicity refinement none|fixed|image (default: image)");
|
||||
logger.Info(" -m<txt> Mosaicity refinement none|image (default: image)");
|
||||
logger.Info(" -A Anomalous mode (don't merge Friedel pairs)");
|
||||
}
|
||||
|
||||
@@ -67,7 +67,7 @@ int main(int argc, char **argv) {
|
||||
bool anomalous_mode = false;
|
||||
std::optional<int> space_group_number;
|
||||
|
||||
enum class MosaicityRefinementMode { None, Fixed, Image };
|
||||
enum class MosaicityRefinementMode { None, Image };
|
||||
MosaicityRefinementMode mosaicity_refinement_mode = MosaicityRefinementMode::Image;
|
||||
|
||||
float d_min_spot_finding = 1.5;
|
||||
@@ -125,8 +125,6 @@ int main(int argc, char **argv) {
|
||||
case 'm':
|
||||
if (strcmp(optarg, "none") == 0)
|
||||
mosaicity_refinement_mode = MosaicityRefinementMode::None;
|
||||
else if (strcmp(optarg, "fixed") == 0)
|
||||
mosaicity_refinement_mode = MosaicityRefinementMode::Fixed;
|
||||
else if (strcmp(optarg, "image") == 0)
|
||||
mosaicity_refinement_mode = MosaicityRefinementMode::Image;
|
||||
else {
|
||||
@@ -433,13 +431,8 @@ int main(int argc, char **argv) {
|
||||
case MosaicityRefinementMode::None:
|
||||
scale_opts.refine_mosaicity = false;
|
||||
break;
|
||||
case MosaicityRefinementMode::Fixed:
|
||||
scale_opts.refine_mosaicity = true;
|
||||
scale_opts.per_image_mosaicity = false;
|
||||
break;
|
||||
case MosaicityRefinementMode::Image:
|
||||
scale_opts.refine_mosaicity = true;
|
||||
scale_opts.per_image_mosaicity = true;
|
||||
break;
|
||||
}
|
||||
if (space_group)
|
||||
@@ -453,10 +446,8 @@ int main(int argc, char **argv) {
|
||||
double scale_time = std::chrono::duration<double>(scale_end - scale_start).count();
|
||||
|
||||
if (scale_result) {
|
||||
logger.Info("Scaling completed in {:.2f} s ({} unique reflections, {} images)",
|
||||
scale_time,
|
||||
scale_result->merged.size(),
|
||||
scale_result->image_ids.size());
|
||||
logger.Info("Scaling completed in {:.2f} s ({} unique reflections)",
|
||||
scale_time, scale_result->merged.size());
|
||||
|
||||
// Print resolution-shell statistics table
|
||||
{
|
||||
@@ -499,13 +490,10 @@ int main(int argc, char **argv) {
|
||||
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_g[i] << "\n";
|
||||
for (size_t i = 0; i < scale_result->mosaicity_deg.size(); ++i) {
|
||||
img_file << i << " " << scale_result->mosaicity_deg[i] << " " << scale_result->image_scale_g[i] << "\n";
|
||||
}
|
||||
img_file.close();
|
||||
logger.Info("Wrote {} image records to {}", scale_result->image_ids.size(), img_path);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Reference in New Issue
Block a user