diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 71d9689c..0bdcbcdf 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -117,36 +117,36 @@ struct IntensityResidual { IntensityResidual(const Reflection& r, double sigma_obs, double wedge_deg, bool refine_partiality) : Iobs_(static_cast(r.I)), inv_sigma_(SafeInv(sigma_obs, 1.0)), - delta_phi_rad_(static_cast(r.delta_phi) * M_PI / 180.0), + delta_phi_(r.delta_phi), lp_(SafeInv(static_cast(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 - 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& ob const int nimg = static_cast(imgIdToSlot.size()); const int nhkl = static_cast(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 k(nimg, 1.0); + std::vector g(nimg, 1.0); std::vector Itrue(nhkl, 0.0); - auto deg2rad = [](double deg) { return deg * M_PI / 180.0; }; - // Mosaicity: always per-image - std::vector mosaicity_rad(nimg, deg2rad(opt.mosaicity_init_deg)); + std::vector mosaicity(nimg, opt.mosaicity_init_deg); // Initialize Itrue from per-HKL median of observed intensities { @@ -286,8 +284,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& 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& ob for (int i = 0; i < nimg; ++i) { auto* rcost = new ceres::AutoDiffCostFunction( 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& 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(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(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(o.r->I)) / o.sigma; - sum_r2 += resid * resid; - } - - const double gof2 = (n_obs > n_params) ? sum_r2 / static_cast(n_obs - n_params) : 1.0; - out.gof2 = gof2; - - // --- Covariance: σ(I_true) from (J^T W J)^{-1} scaled by GoF² --- std::vector slotToHKL(nhkl); for (const auto& kv : hklToSlot) { slotToHKL[kv.second] = kv.first; @@ -399,7 +348,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& 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::quiet_NaN(); + out.merged[h].sigma = 0.0; out.merged[h].d = 0.0; } @@ -420,26 +369,5 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob } } - ceres::Covariance::Options cov_options; - cov_options.algorithm_type = ceres::SPARSE_QR; - - ceres::Covariance covariance(cov_options); - - std::vector> 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; } diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 1e0271fd..65f5f310 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -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 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 merged; - std::vector image_scale_k; + std::vector image_scale_g; std::vector image_ids; // One mosaicity value per image (degrees). diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 40a15f3d..c3c359a0 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -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);