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

ScaleMerge: Simplify by removing log residual

Browse Source
This commit is contained in:
leonarski_f
2026-02-20 09:28:50 +01:00
parent eefea7f36b
commit 43c3402198
3 changed files with 18 additions and 103 deletions
Download Patch File Download Diff File Expand all files Collapse all files

95
image_analysis/scale_merge/ScaleAndMerge.cpp
View File

@@ -115,59 +115,6 @@ namespace {
return key;
}
/// CrystFEL-like log-scaling residual
///
/// residual = w * [ ln(I_obs) - ln(G) - ln(partiality) - ln(lp) - ln(I_true) ]
///
/// Only observations with I_obs > 0 should be fed in (the caller skips the rest).
/// G and I_true are constrained to be positive via Ceres lower bounds.
struct LogIntensityResidual {
LogIntensityResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality)
: log_Iobs_(std::log(std::max(static_cast<double>(r.I), 1e-30))),
weight_(SafeInv(sigma_obs / std::max(static_cast<double>(r.I), 1e-30), 1.0)),
delta_phi_(r.delta_phi_deg),
log_lp_(std::log(std::max(SafeInv(static_cast<double>(r.rlp), 1.0), 1e-30))),
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 G,
const T *const mosaicity,
const T *const Itrue,
T *residual) const {
T partiality;
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_);
}
// Clamp partiality away from zero so log is safe
const T min_p = T(1e-30);
if (partiality < min_p)
partiality = min_p;
// ln(I_pred) = ln(G) + ln(partiality) + ln(lp) + ln(Itrue)
const T log_Ipred = ceres::log(G[0]) + ceres::log(partiality) + T(log_lp_) + ceres::log(Itrue[0]);
residual[0] = (log_Ipred - T(log_Iobs_)) * T(weight_);
return true;
}
double log_Iobs_;
double weight_; // w_i ≈ I_obs / sigma_obs (relative weight in log-space)
double delta_phi_;
double log_lp_;
double half_wedge_;
double c1_;
double partiality_;
bool refine_partiality_;
};
struct IntensityResidual {
IntensityResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality)
: Iobs_(static_cast<double>(r.I)),
@@ -362,32 +309,15 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &ob
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));
if (opt.log_scaling_residual) {
// Log residual requires positive I_obs
if (o.r->I <= 0.0f)
continue;
auto *cost = new ceres::AutoDiffCostFunction<LogIntensityResidual, 1, 1, 1, 1>(
new LogIntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality));
problem.AddResidualBlock(cost,
nullptr,
&g[o.img_slot],
&mosaicity[mos_slot],
&Itrue[o.hkl_slot]);
is_valid_hkl_slot[o.hkl_slot] = true;
} else {
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],
&Itrue[o.hkl_slot]);
is_valid_hkl_slot[o.hkl_slot] = true;
}
problem.AddResidualBlock(cost,
nullptr,
&g[o.img_slot],
&mosaicity[mos_slot],
&Itrue[o.hkl_slot]);
is_valid_hkl_slot[o.hkl_slot] = true;
}
if (opt.smoothen_g) {
@@ -418,15 +348,6 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &ob
for (int i = 0; i < nimg; ++i)
problem.SetParameterLowerBound(&g[i], 0, 1e-12);
// For log residual, Itrue must stay positive
if (opt.log_scaling_residual) {
for (int h = 0; h < nhkl; ++h) {
if (is_valid_hkl_slot[h])
problem.SetParameterLowerBound(&Itrue[h], 0, 1e-12);
}
}
// Mosaicity refinement + bounds
if (!opt.refine_mosaicity) {
for (int i = 0; i < (opt.per_image_mosaicity ? nimg : 1); ++i)

1
image_analysis/scale_merge/ScaleAndMerge.h
View File

@@ -42,7 +42,6 @@ struct ScaleMergeOptions {
// --- Optional: regularize per-image scale k towards 1 (Kabsch-like) ---
bool regularize_scale_to_one = true;
double scale_regularization_sigma = 0.05;
bool log_scaling_residual = false;
double min_partiality_for_merge = 0.2;
bool smoothen_g = true;

25
tools/jfjoch_process.cpp
View File

@@ -42,7 +42,6 @@ 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(" -L Use log-scaling residual");
logger.Info(" -m<txt> Mosaicity refinement none|fixed|image (default: image)");
logger.Info(" -A Anomalous mode (don't merge Friedel pairs)");
}
@@ -68,13 +67,11 @@ int main(int argc, char **argv) {
bool anomalous_mode = false;
std::optional<int> space_group_number;
double resolution_limit = 0.0;
bool log_residual = false;
enum class MosaicityRefinementMode { None, Fixed, Image };
MosaicityRefinementMode mosaicity_refinement_mode = MosaicityRefinementMode::Image;
float d_high = 1.5;
float d_min_spot_finding = 1.5;
std::optional<float> d_min_scale_merge;
if (argc == 1) {
print_usage(logger);
@@ -82,7 +79,7 @@ int main(int argc, char **argv) {
}
int opt;
while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MLm:AD:")) != -1) {
while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:Mm:AD:")) != -1) {
switch (opt) {
case 'o':
output_prefix = optarg;
@@ -100,7 +97,7 @@ int main(int argc, char **argv) {
verbose = true;
break;
case 'd':
resolution_limit = atof(optarg);
d_min_spot_finding = atof(optarg);
break;
case 'R':
rotation_indexing = true;
@@ -113,8 +110,8 @@ int main(int argc, char **argv) {
refine_beam_center = false;
break;
case 'D':
d_high = atof(optarg);
logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_high);
d_min_scale_merge = atof(optarg);
logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_spot_finding);
break;
case 'S':
space_group_number = atoi(optarg);
@@ -122,9 +119,6 @@ int main(int argc, char **argv) {
case 'M':
run_scaling = true;
break;
case 'L':
log_residual = true;
break;
case 'A':
anomalous_mode = true;
break;
@@ -213,7 +207,9 @@ int main(int argc, char **argv) {
SpotFindingSettings spot_settings;
spot_settings.enable = true;
spot_settings.indexing = true;
spot_settings.high_resolution_limit = d_high;
spot_settings.high_resolution_limit = d_min_spot_finding;
if (d_min_scale_merge > 0)
spot_settings.high_resolution_limit = d_min_spot_finding;
// Initialize Analysis Components
PixelMask pixel_mask = dataset->pixel_mask;
@@ -430,9 +426,8 @@ int main(int argc, char **argv) {
scale_opts.refine_mosaicity = true;
scale_opts.max_num_iterations = 500;
scale_opts.max_solver_time_s = 240.0; // generous cutoff for now
scale_opts.log_scaling_residual = log_residual;
scale_opts.merge_friedel = !anomalous_mode;
scale_opts.d_min_limit_A = d_high;
scale_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0);
switch (mosaicity_refinement_mode) {
case MosaicityRefinementMode::None: