mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests 1 Actions Packages Projects Releases 54 Wiki Activity

ScaleAndMerge: Fix symbol names k --> g, remove setting first image scale factor to 1 (it is too weak), remove radians

Browse Source
This commit is contained in:
leonarski_f
2026-02-09 18:00:23 +01:00
parent 684c3df204
commit e0fd62df2b
3 changed files with 28 additions and 102 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/ScaleAndMerge.cpp
+25 -97
View File
@@ -117,36 +117,36 @@ struct IntensityResidual {
IntensityResidual(const Reflection& r, double sigma_obs, double wedge_deg, bool refine_partiality)
: Iobs_(static_cast<double>(r.I)),
inv_sigma_(SafeInv(sigma_obs, 1.0)),
delta_phi_rad_(static_cast<double>(r.delta_phi) * M_PI / 180.0),
delta_phi_(r.delta_phi),
lp_(SafeInv(static_cast<double>(r.rlp), 1.0)),
half_wedge_rad_(wedge_deg * M_PI / 180.0 / 2.0),
half_wedge_(wedge_deg / 2.0),
c1_(r.zeta / std::sqrt(2.0)),
partiality_(r.partiality),
refine_partiality_(refine_partiality) {}
template<typename T>
bool operator()(const T* const k,
const T* const mosaicity_rad,
bool operator()(const T* const G,
const T* const mosaicity,
const T* const Itrue,
T* residual) const {
T partiality;
if (refine_partiality_ && mosaicity_rad[0] != 0.0) {
const T arg_plus = T(delta_phi_rad_ + half_wedge_rad_) * T(c1_) / mosaicity_rad[0];
const T arg_minus = T(delta_phi_rad_ - half_wedge_rad_) * T(c1_) / mosaicity_rad[0];
if (refine_partiality_ && mosaicity[0] != 0.0) {
const T arg_plus = T(delta_phi_ + half_wedge_) * T(c1_) / mosaicity[0];
const T arg_minus = T(delta_phi_ - half_wedge_) * T(c1_) / mosaicity[0];
partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0);
} else
partiality = T(partiality_);
const T Ipred = k[0] * partiality * T(lp_) * Itrue[0];
const T Ipred = G[0] * partiality * T(lp_) * Itrue[0];
residual[0] = (Ipred - T(Iobs_)) * T(inv_sigma_);
return true;
}
double Iobs_;
double inv_sigma_;
double delta_phi_rad_;
double delta_phi_;
double lp_;
double half_wedge_rad_;
double half_wedge_;
double c1_;
double partiality_;
bool refine_partiality_;
@@ -241,20 +241,18 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& ob
const int nimg = static_cast<int>(imgIdToSlot.size());
const int nhkl = static_cast<int>(hklToSlot.size());
out.image_scale_k.assign(nimg, 1.0);
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> k(nimg, 1.0);
std::vector<double> g(nimg, 1.0);
std::vector<double> Itrue(nhkl, 0.0);
auto deg2rad = [](double deg) { return deg * M_PI / 180.0; };
// Mosaicity: always per-image
std::vector<double> mosaicity_rad(nimg, deg2rad(opt.mosaicity_init_deg));
std::vector<double> mosaicity(nimg, opt.mosaicity_init_deg);
// Initialize Itrue from per-HKL median of observed intensities
{
@@ -286,8 +284,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& ob
problem.AddResidualBlock(cost,
nullptr, // no loss function
&k[o.img_slot],
&mosaicity_rad[o.img_slot],
&g[o.img_slot],
&mosaicity[o.img_slot],
&Itrue[o.hkl_slot]);
}
@@ -296,44 +294,25 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& ob
for (int i = 0; i < nimg; ++i) {
auto* rcost = new ceres::AutoDiffCostFunction<ScaleRegularizationResidual, 1, 1>(
new ScaleRegularizationResidual(opt.scale_regularization_sigma));
problem.AddResidualBlock(rcost, nullptr, &k[i]);
problem.AddResidualBlock(rcost, nullptr, &g[i]);
}
}
// Fix gauge freedom
if (opt.fix_first_image_scale && nimg > 0) {
k[0] = 1.0;
problem.SetParameterBlockConstant(&k[0]);
}
// Mosaicity refinement + bounds
if (!opt.refine_mosaicity) {
for (int i = 0; i < nimg; ++i)
problem.SetParameterBlockConstant(&mosaicity_rad[i]);
problem.SetParameterBlockConstant(&mosaicity[i]);
} else {
const double min_rad = deg2rad(opt.mosaicity_min_deg);
const double max_rad = deg2rad(opt.mosaicity_max_deg);
for (int i = 0; i < nimg; ++i) {
problem.SetParameterLowerBound(&mosaicity_rad[i], 0, min_rad);
problem.SetParameterUpperBound(&mosaicity_rad[i], 0, max_rad);
}
}
// Force k[i] > 0 for all images (scale factors must be positive)
{
constexpr double k_floor = 1e-8;
for (int i = 0; i < nimg; ++i) {
if (!(opt.fix_first_image_scale && i == 0)) {
problem.SetParameterLowerBound(&k[i], 0, k_floor);
}
problem.SetParameterLowerBound(&mosaicity[i], 0, opt.mosaicity_min_deg);
problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg);
}
}
ceres::Solver::Options options;
options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
options.minimizer_progress_to_stdout = false;
options.logging_type = ceres::LoggingType::SILENT;
options.minimizer_progress_to_stdout = true;
options.max_num_iterations = opt.max_num_iterations;
options.max_solver_time_in_seconds = opt.max_solver_time_s;
@@ -342,52 +321,22 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& ob
// --- Export per-image results ---
for (int i = 0; i < nimg; ++i)
out.image_scale_k[i] = k[i];
out.image_scale_g[i] = g[i];
out.mosaicity_deg.resize(nimg);
for (int i = 0; i < nimg; ++i)
out.mosaicity_deg[i] = mosaicity_rad[i] * 180.0 / M_PI;
out.mosaicity_deg[i] = mosaicity[i];
// --- Compute goodness-of-fit (reduced chi-squared) ---
const int n_obs = static_cast<int>(obs.size());
// Count free parameters: nhkl Itrue + per-image (k + mosaicity) minus fixed ones
int n_params = nhkl;
for (int i = 0; i < nimg; ++i) {
if (!(opt.fix_first_image_scale && i == 0))
n_params += 1; // k
n_params += 1; // k
if (opt.refine_mosaicity)
n_params += 1; // mosaicity
}
const double half_wedge_rad = opt.wedge_deg * M_PI / 180.0 / 2.0;
double sum_r2 = 0.0;
for (const auto& o : obs) {
const int i = o.img_slot;
const int h = o.hkl_slot;
const double ki = k[i];
const double mosa_i = mosaicity_rad[i];
const double c2 = o.r->zeta / (std::sqrt(2.0) * mosa_i);
const double delta_phi_rad = static_cast<double>(o.r->delta_phi) * M_PI / 180.0;
const double partiality = refine_partiality
? (std::erf((delta_phi_rad + half_wedge_rad) * c2)
- std::erf((delta_phi_rad - half_wedge_rad) * c2)) / 2.0
: o.r->partiality;
const double lp = SafeInv(o.r->rlp, 1.0);
const double Ipred = ki * partiality * lp * Itrue[h];
const double resid = (Ipred - static_cast<double>(o.r->I)) / o.sigma;
sum_r2 += resid * resid;
}
const double gof2 = (n_obs > n_params) ? sum_r2 / static_cast<double>(n_obs - n_params) : 1.0;
out.gof2 = gof2;
// --- Covariance: σ(I_true) from (J^T W J)^{-1} scaled by GoF² ---
std::vector<HKLKey> slotToHKL(nhkl);
for (const auto& kv : hklToSlot) {
slotToHKL[kv.second] = kv.first;
@@ -399,7 +348,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& ob
out.merged[h].k = slotToHKL[h].k;
out.merged[h].l = slotToHKL[h].l;
out.merged[h].I = Itrue[h];
out.merged[h].sigma = std::numeric_limits<double>::quiet_NaN();
out.merged[h].sigma = 0.0;
out.merged[h].d = 0.0;
}
@@ -420,26 +369,5 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection>& ob
}
}
ceres::Covariance::Options cov_options;
cov_options.algorithm_type = ceres::SPARSE_QR;
ceres::Covariance covariance(cov_options);
std::vector<std::pair<const double*, const double*>> covariance_blocks;
covariance_blocks.reserve(nhkl);
for (int h = 0; h < nhkl; ++h) {
covariance_blocks.emplace_back(&Itrue[h], &Itrue[h]);
}
if (covariance.Compute(covariance_blocks, &problem)) {
for (int h = 0; h < nhkl; ++h) {
double var_I = 0.0;
covariance.GetCovarianceBlock(&Itrue[h], &Itrue[h], &var_I);
// σ(I) = I * sqrt( var(I) * GoF² )
out.merged[h].sigma = std::sqrt(var_I * gof2);
}
}
return out;
}
image_analysis/scale_merge/ScaleAndMerge.h
+2 -4
View File
@@ -17,8 +17,6 @@ struct ScaleMergeOptions {
double image_number_rounding = 1.0;
double min_sigma = 1e-3;
bool fix_first_image_scale = true;
// Symmetry canonicalization of HKL prior to merging/scaling.
// If not set, the routine uses raw HKL as-is.
std::optional<gemmi::SpaceGroup> space_group;
@@ -40,7 +38,7 @@ struct ScaleMergeOptions {
double mosaicity_max_deg = 2.0;
// --- Optional: regularize per-image scale k towards 1 (Kabsch-like) ---
bool regularize_scale_to_one = false;
bool regularize_scale_to_one = true;
double scale_regularization_sigma = 0.05;
};
@@ -55,7 +53,7 @@ struct MergedReflection {
struct ScaleMergeResult {
std::vector<MergedReflection> merged;
std::vector<double> image_scale_k;
std::vector<double> image_scale_g;
std::vector<int> image_ids;
// One mosaicity value per image (degrees).
tools/jfjoch_process.cpp
+1 -1
View File
@@ -439,7 +439,7 @@ int main(int argc, char **argv) {
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";
<< scale_result->image_scale_g[i] << "\n";
}
img_file.close();
logger.Info("Wrote {} image records to {}", scale_result->image_ids.size(), img_path);