From 102a2a7c814c8a54a14c95d888e3ed4220ef7095 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Mon, 15 Jun 2026 13:34:54 +0200 Subject: [PATCH] PixelRefine: strip experimental env knobs (orientation, sweep, Lorentz, ML census) 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 --- image_analysis/IndexAndRefine.cpp | 4 +- .../pixel_refinement/PixelRefine.cpp | 405 +++--------------- image_analysis/pixel_refinement/PixelRefine.h | 29 +- 3 files changed, 62 insertions(+), 376 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index a37a6962..c4c0e831 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -448,9 +448,7 @@ bool IndexAndRefine::PixelRefineIntegrate(DataMessage &msg, // Signal-box radius from the shared integration setting (same knob as BraggIntegrate2D). prd.shoebox_radius = static_cast(std::lround(experiment.GetBraggIntegrationSettings().GetR1())); - if (const char *m = std::getenv("PR_RMULT")) prd.r1_multiplier = std::stod(m); // TEMP: Term-2 R1 multiplier sweep - if (std::getenv("PR_ORIENT")) prd.refine_orientation = true; // A/B: per-image orientation refinement - if (std::getenv("PR_SWEEP")) prd.sweep_orientation = true; // A/B: orientation + cell-scale sweep + if (const char *m = std::getenv("PR_RMULT")) prd.r1_multiplier = std::stod(m); // Term-2 R1 multiplier (default 6) std::vector buffer; const uint8_t *ptr = image.GetUncompressedPtr(buffer); diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index b2f2161f..4f45c7aa 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -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 - 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 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 e_obs_recip; - ObservedRecip(beam, distance_mm, detector_rot, obs.x, obs.y, pixel_size, inv_lambda, e_obs_recip); - const Eigen::Matrix 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 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 - 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 -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::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(std::lround(px)); - const int cy = static_cast(std::lround(py)); - const int outer = radius + 1; - if (cx - outer < 0 || cy - outer < 0 || - cx + outer >= static_cast(xpixel) || cy + outer >= static_cast(ypixel)) - return qnan; - double sig = 0.0; - int nsig = 0; - std::vector 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(xpixel) * y + x]; - if (raw == std::numeric_limits::max()) - return qnan; - if (std::is_signed_v && raw == std::numeric_limits::min()) - return qnan; - const double v = static_cast(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(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(data.ewald_dist_cutoff), - .max_hkl = 100, - .centering = data.centering, - .bandwidth_sigma = static_cast(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; - double r_max = 0.0, r_min = std::numeric_limits::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(m.h) + B * static_cast(m.k) - + C * static_cast(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(std::lround(data.sweep_max_deg * M_PI / 180.0 * r_max)), 1, 25); - const int n_scale = std::clamp( - static_cast(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(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 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 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::max()) - continue; - if (std::is_signed_v && (image[npixel] == std::numeric_limits::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(x), static_cast(y), - static_cast(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::max()) + continue; + if (std::is_signed_v && (image[npixel] == std::numeric_limits::min())) + continue; + g.pixels.push_back({static_cast(x), static_cast(y), + static_cast(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( - 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(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(npix)) / sigma_rad; - oprob.AddResidualBlock(new ceres::AutoDiffCostFunction( - new OrientationRegularizer(w, orient_prior)), nullptr, latt_vec0); - } - if (data.scale_reg_sigma > 0.0) { - const double w = std::sqrt(static_cast(groups.size()) / data.scale_reg_sigma); - oprob.AddResidualBlock(new ceres::AutoDiffCostFunction( - 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*) from the intensity-weighted tangential second moment of // the strong spots, binned by resolution (low res small spots, high res larger). diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index 40dba239..366d8ca3 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -68,28 +68,12 @@ struct PixelRefineData { bool refine_B = true; // refine the per-image B-factor along with G // Term 2 measures the physical tangential width R1, but the *integration* profile must - // be generous (XDS-style: integrate over ~6-10 sigma), or a tight template centred on - // the prediction sits off the ~0.4 px centroid-floor scatter and underestimates the - // intensity. The measured R1 is multiplied by this before use. 1.0 = the raw physical - // width (too tight - regressed R-free on the serial jet); ~6 recovers it (closes the - // CCref gap to box-sum). Tune the final value against R-free. + // be generous (XDS-style: integrate over ~6 sigma) or a tight template centred on the + // prediction sits off the ~0.4 px centroid-floor scatter and underestimates the + // intensity (validated on the jet R-free: 0.34 with the raw width -> 0.26 at x6). The + // measured R1 is multiplied by this before use. double r1_multiplier = 6.0; - // Optional per-image orientation refinement (env-gated, off by default): the - // pre-factored per-pixel step that refines the crystal orientation against the - // shoebox pixels, regularised to the spot-centroid orientation, before integration. - // Reinstated to A/B its effect on serial-data R-free vs the fixed-geometry default. - bool refine_orientation = false; - double orient_reg_sigma_deg = 1.0; // orientation-prior strength (deg) - double fit_signal_sigma_pix = 1.5; // signal-weighting sigma for the orientation fit (px) - - // Optional reference-driven orientation + per-image cell-scale sweep before prediction - // (env-gated, off by default). Coarse global moves the gradient refinement cannot make, - // plus a radial (cell-scale) DOF the gradient path lacks. See SweepOrientationCell. - bool sweep_orientation = false; - double sweep_max_deg = 0.15; - double sweep_max_cell_frac = 0.003; - // Relative X-ray bandwidth (sigma of dlambda/lambda), e.g. ~0.004 for a 1% FWHM // DMM, ~1e-4 for Si(111). Adds a resolution-dependent radial broadening to R[0]. // 0 = monochromatic (the term switches off entirely). @@ -136,11 +120,6 @@ class PixelRefine { double beam[2], double &dist_mm, double detector_rot[2], double latt_vec0[3], double latt_vec1[3], double latt_vec2[3]) const; - - // Reference-driven global orientation + cell-scale sweep (see PixelRefineData::sweep_*). - template - void SweepOrientationCell(const T *image, BraggPrediction &prediction, - PixelRefineData &data) const; public: PixelRefine(const DiffractionExperiment &experiment, const std::vector &reference);