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

PixelRefine: strip experimental env knobs (orientation, sweep, Lorentz, ML census)

Browse Source 
Remove the env-gated experiments that A/B'd to dead-ends or are no longer needed,
returning PixelRefine to the clean factored Terms 1+2 plus the one validated keeper
(r1_multiplier, default 6):
- PR_ORIENT (per-image orientation refinement): R-free no-op (0.2618 vs 0.2625) -
  XtalOptimizer's orientation is already optimal. Removes ShoeboxResidual,
  OrientationRegularizer, PixelObs::weight, the refinement block and its fields.
- PR_SWEEP (orientation + cell-scale sweep): R-free no-op, degraded high-res CC1/2
  (per-image overfit). Removes SweepOrientationCell and its fields.
- PR_LORENTZ (rotation Lorentz/zeta): hurt both directions (the factored partiality
  already subsumes it); was already reverted.
- PR_MLCENSUS (multi-lattice census in AnalyzeIndexing): served its purpose (~3-5%
  of jet frames are multi-lattice; shelved).

PR_RMULT (the validated Term-2 multiplier knob) is kept. Defaults unchanged:
crystal 2 / jet / hybrid -R -r pixelrefine all reproduce.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-06-15 13:34:54 +02:00
co-authored by Claude Opus 4.8
parent 6f733d74c2
commit 102a2a7c81
3 changed files with 62 additions and 376 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/pixel_refinement/PixelRefine.cpp
+57 -348
View File
@@ -24,7 +24,6 @@ struct PixelObs {
double x, y; // detector pixel coordinate
double Iobs; // raw pixel value (signal + background)
double Ibkg; // local background estimate (per-shoebox level, raw counts)
double weight; // 1/sigma_pixel (used only by the optional orientation refinement)
};
// One reflection together with the pixels of its shoebox.
@@ -337,69 +336,6 @@ struct IntensityResidual {
double J, inv_sigma, partiality, pol, I_ref, inv_4d2;
};
// Per-shoebox per-pixel forward-model cost (raw counts), used ONLY by the optional
// orientation refinement: I_pred = G*Itrue*B_term*P_radial*P_tang*pol + Ibkg. The
// expensive node geometry is computed once per reflection; the cheap ObservedRecip +
// Gaussian profile run per pixel. Shares ObservedRecip/PredictedNode with GeometryProbe.
struct ShoeboxResidual {
ShoeboxResidual(const ReflGroup &g, double lambda, double pixel_size,
gemmi::CrystalSystem symmetry)
: pixels(g.pixels), Itrue(g.Itrue), R_bw_sq(g.R_bw_sq), pol(g.pol),
exp_h(g.h), exp_k(g.k), exp_l(g.l),
inv_lambda(1.0 / lambda), pixel_size(pixel_size), symmetry(symmetry) {}
template<typename T>
bool operator()(const T *const *params, T *residual) const {
// 0 beam 1 dist 2 detector_rot 3 p0 4 p1 5 p2 6 scale 7 B 8 R
const T *beam = params[0]; const T *distance_mm = params[1]; const T *detector_rot = params[2];
const T *p0 = params[3]; const T *p1 = params[4]; const T *p2 = params[5];
const T *scale_factor = params[6]; const T *B = params[7]; const T *R = params[8];
if (R[0] < T(1e-10) || R[1] < T(1e-10))
return false;
Eigen::Matrix<T, 3, 1> e_pred_recip, n_radial;
T q_sq;
if (!PredictedNode(p0, p1, p2, exp_h, exp_k, exp_l, symmetry, inv_lambda,
e_pred_recip, n_radial, q_sq))
return false;
const T B_term = ceres::exp(-B[0] * q_sq / T(4.0));
const T R0_eff_sq = R[0] * R[0] + T(R_bw_sq);
for (size_t i = 0; i < pixels.size(); ++i) {
const PixelObs &obs = pixels[i];
Eigen::Matrix<T, 3, 1> e_obs_recip;
ObservedRecip(beam, distance_mm, detector_rot, obs.x, obs.y, pixel_size, inv_lambda, e_obs_recip);
const Eigen::Matrix<T, 3, 1> delta_q = e_obs_recip - e_pred_recip;
const T eps_radial = delta_q.dot(n_radial);
const T eps_tang_sq = (delta_q - eps_radial * n_radial).squaredNorm();
const T P_radial = ceres::exp(-eps_radial * eps_radial / R0_eff_sq);
const T P_tang = ceres::exp(-eps_tang_sq / (R[1] * R[1])) / (T(M_PI) * R[1] * R[1]);
const T signal = scale_factor[0] * T(Itrue) * B_term * P_radial * P_tang * T(pol);
residual[i] = (signal + T(obs.Ibkg) - T(obs.Iobs)) * T(obs.weight);
}
return true;
}
std::vector<PixelObs> pixels;
const double Itrue, R_bw_sq, pol;
const double exp_h, exp_k, exp_l;
const double inv_lambda, pixel_size;
gemmi::CrystalSystem symmetry;
};
// Anchors the orientation (angle-axis vector) to its prior with a data-scaled weight, so
// the per-image fit can only make a small, signal-supported sub-spot correction.
struct OrientationRegularizer {
OrientationRegularizer(double weight, const double prior[3]) : weight(weight) {
for (int i = 0; i < 3; ++i)
prior_[i] = prior[i];
}
template<typename T>
bool operator()(const T *p0, T *residual) const {
for (int i = 0; i < 3; ++i)
residual[i] = T(weight) * (p0[i] - T(prior_[i]));
return true;
}
double weight;
double prior_[3];
};
} // namespace
PixelRefine::PixelRefine(const DiffractionExperiment &experiment,
@@ -457,144 +393,6 @@ void PixelRefine::BuildParameterBlocks(const PixelRefineData &data,
}
}
// Optional pre-pass (env-gated): a small GLOBAL orientation + uniform cell-scale sweep that
// maximises CC of the box-summed intensities against the reference. Unlike the per-pixel
// orientation refinement it also adjusts a per-image cell scale (a radial degree of freedom),
// and makes coarse global moves the local gradient cannot. Coordinate descent over the three
// Rodrigues axes + cell scale within geometry-derived pixel bounds (highest-resolution spot
// moves ~1 px/step, lowest barely moves). Writes the best orientation/cell into data.latt.
template<class T>
void PixelRefine::SweepOrientationCell(const T *image, BraggPrediction &prediction,
PixelRefineData &data) const {
const int radius = data.shoebox_radius;
const double beam_x = data.geom.GetBeamX_pxl();
const double beam_y = data.geom.GetBeamY_pxl();
const auto qnan = std::numeric_limits<double>::quiet_NaN();
// Box-sum minus local (perimeter) background MEAN, raw counts; NaN off-detector/masked.
auto integrate = [&](double px, double py) -> double {
const int cx = static_cast<int>(std::lround(px));
const int cy = static_cast<int>(std::lround(py));
const int outer = radius + 1;
if (cx - outer < 0 || cy - outer < 0 ||
cx + outer >= static_cast<int>(xpixel) || cy + outer >= static_cast<int>(ypixel))
return qnan;
double sig = 0.0;
int nsig = 0;
std::vector<double> ring;
ring.reserve((2 * outer + 1) * (2 * outer + 1));
for (int y = cy - outer; y <= cy + outer; ++y) {
for (int x = cx - outer; x <= cx + outer; ++x) {
const T raw = image[static_cast<size_t>(xpixel) * y + x];
if (raw == std::numeric_limits<T>::max())
return qnan;
if (std::is_signed_v<T> && raw == std::numeric_limits<T>::min())
return qnan;
const double v = static_cast<double>(raw);
if (std::abs(x - cx) <= radius && std::abs(y - cy) <= radius) {
sig += v; ++nsig;
} else {
ring.push_back(v);
}
}
}
if (ring.size() < 5)
return qnan;
double rsum = 0.0;
for (const double v : ring)
rsum += v;
return sig - nsig * (rsum / static_cast<double>(ring.size()));
};
DiffractionExperiment exp_iter = experiment;
exp_iter.BeamX_pxl(beam_x).BeamY_pxl(beam_y)
.DetectorDistance_mm(data.geom.GetDetectorDistance_mm())
.PoniRot1_rad(data.geom.GetPoniRot1_rad())
.PoniRot2_rad(data.geom.GetPoniRot2_rad());
const BraggPredictionSettings settings{
.high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
.ewald_dist_cutoff = static_cast<float>(data.ewald_dist_cutoff),
.max_hkl = 100,
.centering = data.centering,
.bandwidth_sigma = static_cast<float>(data.bandwidth)
};
const int nrefl = prediction.Calc(exp_iter, data.latt, settings);
const auto &predicted = prediction.GetReflections();
struct Matched { int h, k, l; double refI; };
std::vector<Matched> matched;
double r_max = 0.0, r_min = std::numeric_limits<double>::max();
for (int i = 0; i < nrefl; ++i) {
const auto &r = predicted[i];
const auto it = reference_data.find(hkl_key_generator(r));
if (it == reference_data.end())
continue;
matched.push_back({r.h, r.k, r.l, it->second});
const double dx = r.predicted_x - beam_x, dy = r.predicted_y - beam_y;
const double rad = std::sqrt(dx * dx + dy * dy);
r_max = std::max(r_max, rad);
r_min = std::min(r_min, rad);
}
if (matched.size() < 20 || r_min <= 1.0 || r_max <= r_min)
return; // too little to anchor a meaningful sweep
auto score = [&](const CrystalLattice &L) -> double {
const Coord A = L.Astar(), B = L.Bstar(), C = L.Cstar();
double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
int n = 0;
for (const auto &m : matched) {
const Coord g = A * static_cast<float>(m.h) + B * static_cast<float>(m.k)
+ C * static_cast<float>(m.l);
const auto [x, y] = data.geom.RecipToDetector(g);
if (!std::isfinite(x) || !std::isfinite(y))
continue;
const double I = integrate(x, y);
if (!std::isfinite(I))
continue;
sx += I; sy += m.refI; sxx += I * I; syy += m.refI * m.refI; sxy += I * m.refI; ++n;
}
if (n < 10)
return -2.0;
const double nd = n;
const double cov = sxy - sx * sy / nd, vx = sxx - sx * sx / nd, vy = syy - sy * sy / nd;
return (vx > 0.0 && vy > 0.0) ? cov / std::sqrt(vx * vy) : -2.0;
};
const double step = 1.0 / r_max;
const int n_rot = std::clamp(
static_cast<int>(std::lround(data.sweep_max_deg * M_PI / 180.0 * r_max)), 1, 25);
const int n_scale = std::clamp(
static_cast<int>(std::lround(data.sweep_max_cell_frac * r_max)), 1, 25);
const Coord axes[3] = {Coord(1, 0, 0), Coord(0, 1, 0), Coord(0, 0, 1)};
CrystalLattice best = data.latt;
double best_cc = score(best);
for (int round = 0; round < 2; ++round) {
for (const auto &axis : axes) {
CrystalLattice axis_best = best;
double axis_cc = best_cc;
for (int i = -n_rot; i <= n_rot; ++i) {
if (i == 0) continue;
CrystalLattice cand = best.Multiply(RotMatrix(static_cast<float>(i * step), axis));
const double cc = score(cand);
if (cc > axis_cc) { axis_cc = cc; axis_best = cand; }
}
best = axis_best; best_cc = axis_cc;
}
CrystalLattice scale_best = best;
double scale_cc = best_cc;
for (int i = -n_scale; i <= n_scale; ++i) {
if (i == 0) continue;
const double s = 1.0 / (1.0 + i * step);
CrystalLattice cand = best.Multiply(gemmi::Mat33(s, 0, 0, 0, s, 0, 0, 0, s));
const double cc = score(cand);
if (cc > scale_cc) { scale_cc = cc; scale_best = cand; }
}
best = scale_best; best_cc = scale_cc;
}
data.latt = best;
}
template<class T>
void PixelRefine::Run(const T *image,
BraggPrediction &prediction,
@@ -602,10 +400,6 @@ void PixelRefine::Run(const T *image,
data.solved = false;
data.reflections.clear();
// Optional reference-driven orientation + cell-scale sweep before prediction (env-gated).
if (data.sweep_orientation)
SweepOrientationCell(image, prediction, data);
const double lambda = data.geom.GetWavelength_A();
const double pixel_size = data.geom.GetPixelSize_mm();
@@ -649,154 +443,69 @@ void PixelRefine::Run(const T *image,
double latt_vec0[3], latt_vec1[3], latt_vec2[3];
BuildParameterBlocks(data, beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2);
// ---- 1-2. Predict shoeboxes + collect pixels (a lambda so it can be re-run after
// an orientation refinement moves the predicted positions). ----------
// A spot-core mask over ALL predictions keeps each reflection's background ring from
// picking up a neighbour's signal. Pixels carry a fit weight (background-limited
// variance, signal-weighted toward the predicted centre) used only by the optional
// orientation refinement - the factored integration weights by v = Ibkg directly.
// ---- 1. Predict shoeboxes for the current geometry ------------------------
DiffractionExperiment exp_iter = experiment;
exp_iter.BeamX_pxl(data.geom.GetBeamX_pxl())
.BeamY_pxl(data.geom.GetBeamY_pxl())
.DetectorDistance_mm(data.geom.GetDetectorDistance_mm())
.PoniRot1_rad(data.geom.GetPoniRot1_rad())
.PoniRot2_rad(data.geom.GetPoniRot2_rad());
const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction);
// ---- 2. Collect per-reflection shoebox pixels + local background ----------
// GetReflections() returns the full pre-sized buffer; only the first nrefl
// entries are valid for this image. A spot-core mask over ALL predictions keeps
// each reflection's background ring from picking up a neighbour's signal.
const auto &predicted = prediction.GetReflections();
const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius);
std::vector<ReflGroup> groups;
auto build_groups = [&]() {
groups.clear();
DiffractionExperiment exp_iter = experiment;
exp_iter.BeamX_pxl(data.geom.GetBeamX_pxl())
.BeamY_pxl(data.geom.GetBeamY_pxl())
.DetectorDistance_mm(data.geom.GetDetectorDistance_mm())
.PoniRot1_rad(data.geom.GetPoniRot1_rad())
.PoniRot2_rad(data.geom.GetPoniRot2_rad());
const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction);
const auto &predicted = prediction.GetReflections();
const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius);
for (int ri = 0; ri < nrefl; ++ri) {
const auto &refl = predicted[ri];
const auto hkl = hkl_key_generator(refl);
if (!reference_data.contains(hkl))
continue;
double Ibkg = 0.0;
if (!EstimateLocalBackground(image, spot_mask, xpixel, ypixel,
refl.predicted_x, refl.predicted_y,
radius, bkg_outer_radius, Ibkg))
continue;
ReflGroup g;
g.h = refl.h;
g.k = refl.k;
g.l = refl.l;
g.d = refl.d;
g.Itrue = reference_data[hkl];
g.R_bw_sq = bandwidth_radial_sq(refl.d);
g.pol = polarization(refl.predicted_x, refl.predicted_y);
g.Ibkg = Ibkg;
g.predicted_x = refl.predicted_x;
g.predicted_y = refl.predicted_y;
const auto box = ShoeboxBounds(refl.predicted_x, refl.predicted_y, radius, xpixel, ypixel);
for (int y = box.min_y; y <= box.max_y; ++y) {
for (int x = box.min_x; x <= box.max_x; ++x) {
const size_t npixel = xpixel * y + x;
if (image[npixel] == std::numeric_limits<T>::max())
continue;
if (std::is_signed_v<T> && (image[npixel] == std::numeric_limits<T>::min()))
continue;
double weight = 1.0 / std::sqrt(std::max(Ibkg, 1.0));
if (data.fit_signal_sigma_pix > 0.0) {
const double dx = x - refl.predicted_x, dy = y - refl.predicted_y;
const double s2 = data.fit_signal_sigma_pix * data.fit_signal_sigma_pix;
weight *= std::exp(-0.5 * (dx * dx + dy * dy) / s2);
}
g.pixels.push_back({static_cast<double>(x), static_cast<double>(y),
static_cast<double>(image[npixel]), Ibkg, weight});
}
for (int ri = 0; ri < nrefl; ++ri) {
const auto &refl = predicted[ri];
const auto hkl = hkl_key_generator(refl);
if (!reference_data.contains(hkl))
continue;
// Local flat background from the ring around the shoebox (raw counts). If we
// cannot estimate a clean local background the reflection is dropped, exactly
// as BraggIntegrate2D marks it unobserved when too few background pixels survive.
double Ibkg = 0.0;
if (!EstimateLocalBackground(image, spot_mask, xpixel, ypixel,
refl.predicted_x, refl.predicted_y,
radius, bkg_outer_radius, Ibkg))
continue;
ReflGroup g;
g.h = refl.h;
g.k = refl.k;
g.l = refl.l;
g.d = refl.d;
g.Itrue = reference_data[hkl];
g.R_bw_sq = bandwidth_radial_sq(refl.d);
g.pol = polarization(refl.predicted_x, refl.predicted_y);
g.Ibkg = Ibkg;
g.predicted_x = refl.predicted_x;
g.predicted_y = refl.predicted_y;
const auto box = ShoeboxBounds(refl.predicted_x, refl.predicted_y, radius, xpixel, ypixel);
for (int y = box.min_y; y <= box.max_y; ++y) {
for (int x = box.min_x; x <= box.max_x; ++x) {
const size_t npixel = xpixel * y + x;
// Skip sentinel (masked / saturated) pixels.
if (image[npixel] == std::numeric_limits<T>::max())
continue;
if (std::is_signed_v<T> && (image[npixel] == std::numeric_limits<T>::min()))
continue;
g.pixels.push_back({static_cast<double>(x), static_cast<double>(y),
static_cast<double>(image[npixel]), Ibkg});
}
if (!g.pixels.empty())
groups.push_back(std::move(g));
}
};
build_groups();
if (!g.pixels.empty())
groups.push_back(std::move(g));
}
if (groups.empty())
return;
// ---- Optional per-image orientation refinement (off by default) -----------------
// The pre-factored per-pixel step: refine the orientation (and a nuisance scale)
// against the shoebox pixels via ShoeboxResidual, regularised to the spot-centroid
// orientation, then re-predict. Dropped from the factored model (geometry fixed to
// XtalOptimizer); behind PixelRefineData::refine_orientation to A/B its R-free effect.
if (data.refine_orientation) {
const double orient_prior[3] = {latt_vec0[0], latt_vec0[1], latt_vec0[2]};
ceres::Problem oprob;
size_t npix = 0;
for (const auto &g : groups) {
auto *cost = new ceres::DynamicAutoDiffCostFunction<ShoeboxResidual>(
new ShoeboxResidual(g, lambda, pixel_size, data.crystal_system));
for (int b : {2, 1, 2, 3, 3, 3, 1, 1, 2})
cost->AddParameterBlock(b);
cost->SetNumResiduals(static_cast<int>(g.pixels.size()));
oprob.AddResidualBlock(cost, nullptr, beam, &dist_mm, detector_rot,
latt_vec0, latt_vec1, latt_vec2,
&data.scale_factor, &data.B_factor, data.R);
npix += g.pixels.size();
}
// Refine only orientation (latt_vec0) + the nuisance scale G; everything else fixed.
oprob.SetParameterBlockConstant(beam);
oprob.SetParameterBlockConstant(&dist_mm);
oprob.SetParameterBlockConstant(detector_rot);
oprob.SetParameterBlockConstant(latt_vec1);
oprob.SetParameterBlockConstant(latt_vec2);
oprob.SetParameterBlockConstant(&data.B_factor);
oprob.SetParameterBlockConstant(data.R);
oprob.SetParameterLowerBound(&data.scale_factor, 0, 0.0);
if (data.orient_reg_sigma_deg > 0.0 && npix > 0) {
const double sigma_rad = std::max(data.orient_reg_sigma_deg * M_PI / 180.0, 1e-9);
const double w = std::sqrt(static_cast<double>(npix)) / sigma_rad;
oprob.AddResidualBlock(new ceres::AutoDiffCostFunction<OrientationRegularizer, 3, 3>(
new OrientationRegularizer(w, orient_prior)), nullptr, latt_vec0);
}
if (data.scale_reg_sigma > 0.0) {
const double w = std::sqrt(static_cast<double>(groups.size()) / data.scale_reg_sigma);
oprob.AddResidualBlock(new ceres::AutoDiffCostFunction<ScalarRegularizer, 1, 1>(
new ScalarRegularizer(w, 1.0)), nullptr, &data.scale_factor);
}
ceres::Solver::Options oopt;
oopt.linear_solver_type = ceres::DENSE_QR;
oopt.logging_type = ceres::LoggingType::SILENT;
oopt.minimizer_progress_to_stdout = false;
oopt.max_solver_time_in_seconds = data.max_time_s;
oopt.num_threads = 1;
ceres::Solver::Summary osum;
ceres::Solve(oopt, &oprob, &osum);
// Write the refined orientation back into data.latt (cell held at latt_vec1/2),
// then re-predict + rebuild groups at the new orientation. Scale is reset; Term 1
// re-fits it properly below.
switch (data.crystal_system) {
case gemmi::CrystalSystem::Tetragonal:
latt_vec1[1] = latt_vec1[0];
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, M_PI/2, M_PI/2);
break;
case gemmi::CrystalSystem::Cubic:
latt_vec1[1] = latt_vec1[0];
latt_vec1[2] = latt_vec1[0];
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, M_PI/2, M_PI/2);
break;
case gemmi::CrystalSystem::Hexagonal:
latt_vec1[1] = latt_vec1[0];
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, M_PI/2, 2.0*M_PI/3.0);
break;
case gemmi::CrystalSystem::Monoclinic:
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, latt_vec2[0], M_PI/2);
break;
case gemmi::CrystalSystem::Orthorhombic:
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, M_PI/2, M_PI/2);
break;
default:
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1,
latt_vec2[0], latt_vec2[1], latt_vec2[2]);
break;
}
data.scale_factor = 1.0;
build_groups();
if (groups.empty())
return;
}
// ---- 3. Term 2: per-resolution tangential profile width R1 ----------------
// R1 = sqrt(2*<eps_t^2>) from the intensity-weighted tangential second moment of
// the strong spots, binned by resolution (low res small spots, high res larger).