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

PixelRefine: Work in progress (Claude)

Browse Source
This commit is contained in:
leonarski_f
2026-06-08 13:44:39 +02:00
parent 66a48c4266
commit 6f6098d508
2 changed files with 360 additions and 243 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/pixel_refinement/PixelRefine.cpp
+293 -243
View File
@@ -238,15 +238,27 @@ struct PixelResidual {
return false;
const T B_term = ceres::exp(-B[0] * q_sq / T(4.0));
const T P_radial = ceres::exp(-eps_radial * eps_radial / (R[0] * R[0]));
// Normalized 2D Gaussian profile in the detector-tangential plane,
// weighted by the reciprocal-space area of the pixel so that the sum of
// the profile over the shoebox is ~1 (intensity-conserving).
const T norm_tang = T(A_recip) / (T(M_PI) * R[1] * R[1]);
const T P_tang = norm_tang * ceres::exp(-eps_tang_sq / (R[1] * R[1]));
// Full 3D reciprocal-space spot density modelled as a separable Gaussian,
// normalized so that its integral is 1 (intensity-conserving):
// radial: g_r(e) = exp(-e^2/R0^2) / (sqrt(pi) R0) [1/A^-1]
// tangential: g_t(e) = exp(-|e|^2/R1^2) / (pi R1^2) [1/A^-2]
// The detector pixel captures the fraction g_t * A_recip of the tangential
// profile (A_recip = reciprocal area the pixel subtends; sum over shoebox
// ~ 1). The radial factor is the still-image partiality: how far this
// reflection sits from the Ewald sphere.
//
// Caveat: a still samples the radial direction at a single offset, so the
// sqrt(pi) R0 normalization makes g_r a density (1/A^-1) rather than a
// dimensionless fraction. The leftover dimensional factor is absorbed by
// the free scale G; completing it physically needs an effective radial
// sampling width from the energy bandwidth + beam divergence (TODO).
const T g_radial = ceres::exp(-eps_radial * eps_radial / (R[0] * R[0]))
/ (ceres::sqrt(T(M_PI)) * R[0]);
const T P_tang = T(A_recip) * 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;
const T signal = scale_factor[0] * T(Itrue) * B_term * g_radial * P_tang;
Ipred = signal + T(Ibkg);
return true;
}
@@ -303,24 +315,14 @@ void PixelRefine::Run(const T *image,
data.solved = false;
data.reflections.clear();
const auto geom = data.geom;
const double lambda = geom.GetWavelength_A();
const double pixel_size = geom.GetPixelSize_mm();
const double distance_mm = geom.GetDetectorDistance_mm();
const double lambda = data.geom.GetWavelength_A();
const double pixel_size = data.geom.GetPixelSize_mm();
// Reciprocal-space area subtended by one pixel (small-angle approximation):
// |dq/dpixel| ~ pixel_size / (lambda * distance), so the area is its square.
// Used to normalize the tangential profile. Good enough for a prototype; a
// proper treatment would use the per-reflection 2-theta dependent Jacobian.
const double A_recip = std::pow(pixel_size / (lambda * distance_mm), 2.0);
// Predict reflection positions with the current lattice/geometry.
const BraggPredictionSettings settings_prediction{
.high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
.max_hkl = 100,
.centering = data.centering
};
prediction.Calc(experiment, data.latt, settings_prediction);
const auto azim_result = profile.GetResult();
const auto azim_std = profile.GetStd();
@@ -331,249 +333,295 @@ void PixelRefine::Run(const T *image,
const double angle_rad = data.angle_deg * M_PI / 180.0;
const int radius = data.shoebox_radius;
// ---- Collect per-reflection shoebox pixels --------------------------------
// Exact reciprocal-space area a 1x1 pixel subtends, |dq/dx x dq/dy|, via
// finite differences of the detector->reciprocal map. This is the Jacobian
// between the curved Ewald-sphere sampling and flat reciprocal space, and it
// is exactly the geometric factor that plays the role of the Lorentz factor
// for stills: where the sphere grazes reciprocal space obliquely, a pixel
// covers more reciprocal volume and the captured fraction grows. It tracks
// the refined geometry because it reads the current data.geom each iteration.
auto recip_area = [&](double x, double y) -> double {
const Coord qx = data.geom.DetectorToRecip(x + 0.5, y) - data.geom.DetectorToRecip(x - 0.5, y);
const Coord qy = data.geom.DetectorToRecip(x, y + 0.5) - data.geom.DetectorToRecip(x, y - 0.5);
return (qx % qy).Length();
};
// Mutable experiment whose geometry is re-synced from the refined data.geom
// before each prediction, so shoeboxes track the refined geometry/cell.
DiffractionExperiment exp_iter = experiment;
// State retained after the loop for the final reflection extraction.
std::vector<ReflGroup> groups;
for (const auto &refl : prediction.GetReflections()) {
const auto hkl = hkl_key_generator(refl);
if (!reference_data.contains(hkl))
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.predicted_x = refl.predicted_x;
g.predicted_y = refl.predicted_y;
const int min_y = std::max<int>(refl.predicted_y - radius, 0);
const int max_y = std::min<int>(refl.predicted_y + radius, ypixel - 1);
const int min_x = std::max<int>(refl.predicted_x - radius, 0);
const int max_x = std::min<int>(refl.predicted_x + radius, xpixel - 1);
for (int y = min_y; y <= max_y; ++y) {
for (int x = min_x; x <= max_x; ++x) {
const size_t npixel = xpixel * y + x;
const int azim_bin = pixel_to_bin[npixel];
// Skip pixels not mapped to an azimuthal bin or carrying a
// sentinel (masked / saturated) value. We assume the pixel mask
// is already applied upstream.
if (azim_bin >= azim_bin_count)
continue;
if (image[npixel] == std::numeric_limits<T>::max())
continue;
if (std::is_signed_v<T> && (image[npixel] == std::numeric_limits<T>::min()))
continue;
const double correction = corrections[npixel];
const double Ibkg = azim_result[azim_bin] * 1.0; // already corrected units
const double Ibkg_sigma = azim_std[azim_bin];
const double raw = static_cast<double>(image[npixel]);
const double Iobs = raw * correction;
// Per-pixel variance: Poisson noise of the corrected counts
// (var(c*N) = c^2 * N = c * Iobs) plus the background spread.
double var = correction * std::max(Iobs, 0.0) + Ibkg_sigma * Ibkg_sigma;
if (!(var > 1.0))
var = 1.0;
PixelObs obs{
.x = static_cast<double>(x),
.y = static_cast<double>(y),
.Iobs = Iobs,
.Ibkg = Ibkg,
.weight = 1.0 / std::sqrt(var),
.A_recip = A_recip,
.angle_rad = angle_rad
};
g.pixels.push_back(obs);
}
}
if (!g.pixels.empty())
groups.push_back(std::move(g));
}
if (groups.empty())
return;
// ---- Set up parameter blocks (mirrors XtalOptimizer for the geometry part) -
double beam[2] = {geom.GetBeamX_pxl(), geom.GetBeamY_pxl()};
double dist_mm = distance_mm;
double detector_rot[2] = {geom.GetPoniRot1_rad(), geom.GetPoniRot2_rad()};
double beam[2] = {0, 0};
double dist_mm = data.geom.GetDetectorDistance_mm();
double detector_rot[2] = {0, 0};
double rot_vec[3] = {1.0, 0.0, 0.0};
if (auto axis = geom.GetRotation()) {
rot_vec[0] = axis->GetAxis().x;
rot_vec[1] = axis->GetAxis().y;
rot_vec[2] = axis->GetAxis().z;
}
double latt_vec0[3] = {0, 0, 0}; // orientation (Rodrigues)
double latt_vec1[3] = {0, 0, 0}; // lengths
double latt_vec2[3] = {0, 0, 0}; // angles (rad)
double beta = data.latt.GetUnitCell().beta;
switch (data.crystal_system) {
case gemmi::CrystalSystem::Orthorhombic:
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
break;
case gemmi::CrystalSystem::Tetragonal:
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
latt_vec1[0] = (latt_vec1[0] + latt_vec1[1]) / 2.0;
break;
case gemmi::CrystalSystem::Cubic:
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
latt_vec1[0] = (latt_vec1[0] + latt_vec1[1] + latt_vec1[2]) / 3.0;
break;
case gemmi::CrystalSystem::Hexagonal:
LatticeToRodriguesAndLengths_Hex(data.latt, latt_vec0, latt_vec1);
break;
case gemmi::CrystalSystem::Monoclinic:
LatticeToRodriguesLengthsBeta_Mono(data.latt, latt_vec0, latt_vec1, beta);
latt_vec2[0] = beta;
break;
default: {
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
const auto uc = data.latt.GetUnitCell();
latt_vec2[0] = uc.alpha * M_PI / 180.0;
latt_vec2[1] = uc.beta * M_PI / 180.0;
latt_vec2[2] = uc.gamma * M_PI / 180.0;
break;
const int n_iter = std::max(1, data.max_iterations);
for (int iter = 0; iter < n_iter; ++iter) {
// ---- 1. Re-sync prediction geometry from the (refined) model ----------
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());
prediction.Calc(exp_iter, data.latt, settings_prediction);
// ---- 2. Collect per-reflection shoebox pixels -------------------------
groups.clear();
for (const auto &refl : prediction.GetReflections()) {
const auto hkl = hkl_key_generator(refl);
if (!reference_data.contains(hkl))
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.predicted_x = refl.predicted_x;
g.predicted_y = refl.predicted_y;
const int min_y = std::max<int>(refl.predicted_y - radius, 0);
const int max_y = std::min<int>(refl.predicted_y + radius, ypixel - 1);
const int min_x = std::max<int>(refl.predicted_x - radius, 0);
const int max_x = std::min<int>(refl.predicted_x + radius, xpixel - 1);
for (int y = min_y; y <= max_y; ++y) {
for (int x = min_x; x <= max_x; ++x) {
const size_t npixel = xpixel * y + x;
const int azim_bin = pixel_to_bin[npixel];
// Skip pixels not mapped to an azimuthal bin or carrying a
// sentinel (masked / saturated) value. We assume the pixel
// mask is already applied upstream.
if (azim_bin >= azim_bin_count)
continue;
if (image[npixel] == std::numeric_limits<T>::max())
continue;
if (std::is_signed_v<T> && (image[npixel] == std::numeric_limits<T>::min()))
continue;
const double correction = corrections[npixel];
const double Ibkg = azim_result[azim_bin]; // already in corrected units
const double Ibkg_sigma = azim_std[azim_bin];
const double raw = static_cast<double>(image[npixel]);
const double Iobs = raw * correction;
// Per-pixel variance: Poisson noise of the corrected counts
// (var(c*N) = c^2 * N = c * Iobs) plus the background spread.
double var = correction * std::max(Iobs, 0.0) + Ibkg_sigma * Ibkg_sigma;
if (!(var > 1.0))
var = 1.0;
PixelObs obs{
.x = static_cast<double>(x),
.y = static_cast<double>(y),
.Iobs = Iobs,
.Ibkg = Ibkg,
.weight = 1.0 / std::sqrt(var),
.A_recip = recip_area(x, y),
.angle_rad = angle_rad
};
g.pixels.push_back(obs);
}
}
if (!g.pixels.empty())
groups.push_back(std::move(g));
}
}
ceres::Problem problem;
if (groups.empty())
return;
for (const auto &g : groups) {
for (const auto &obs : g.pixels) {
auto *cost = new ceres::AutoDiffCostFunction<
PixelResidual, 1, 2, 1, 2, 3, 3, 3, 3, 1, 1, 2>(
new PixelResidual(obs, g.Itrue, lambda, pixel_size,
g.h, g.k, g.l, data.crystal_system));
problem.AddResidualBlock(cost, new ceres::HuberLoss(3.0),
beam, &dist_mm, detector_rot, rot_vec,
latt_vec0, latt_vec1, latt_vec2,
&data.scale_factor, &data.B_factor, data.R);
// ---- 3. Set up parameter blocks (geometry part mirrors XtalOptimizer) -
beam[0] = data.geom.GetBeamX_pxl();
beam[1] = data.geom.GetBeamY_pxl();
dist_mm = data.geom.GetDetectorDistance_mm();
detector_rot[0] = data.geom.GetPoniRot1_rad();
detector_rot[1] = data.geom.GetPoniRot2_rad();
rot_vec[0] = 1.0; rot_vec[1] = 0.0; rot_vec[2] = 0.0;
if (auto axis = data.geom.GetRotation()) {
rot_vec[0] = axis->GetAxis().x;
rot_vec[1] = axis->GetAxis().y;
rot_vec[2] = axis->GetAxis().z;
}
}
data.residual_count = problem.NumResidualBlocks();
// ---- Constrain / bound parameter blocks -----------------------------------
if (!data.refine_orientation)
problem.SetParameterBlockConstant(latt_vec0);
if (!data.refine_unit_cell) {
problem.SetParameterBlockConstant(latt_vec1);
problem.SetParameterBlockConstant(latt_vec2);
} else {
for (int i = 0; i < 3; ++i) {
problem.SetParameterLowerBound(latt_vec1, i, 5.0);
problem.SetParameterUpperBound(latt_vec1, i, 1000.0);
}
if (data.crystal_system != gemmi::CrystalSystem::Monoclinic &&
data.crystal_system != gemmi::CrystalSystem::Triclinic)
problem.SetParameterBlockConstant(latt_vec2);
}
if (!data.refine_beam_center)
problem.SetParameterBlockConstant(beam);
if (!data.refine_distance) {
problem.SetParameterBlockConstant(&dist_mm);
} else {
problem.SetParameterLowerBound(&dist_mm, 0, dist_mm * 0.9);
problem.SetParameterUpperBound(&dist_mm, 0, dist_mm * 1.1);
}
if (!data.refine_detector_angles) {
problem.SetParameterBlockConstant(detector_rot);
} else {
const double rng = 3.0 / 180.0 * M_PI;
for (int i = 0; i < 2; ++i) {
problem.SetParameterLowerBound(detector_rot, i, detector_rot[i] - rng);
problem.SetParameterUpperBound(detector_rot, i, detector_rot[i] + rng);
}
}
if (!data.refine_rotation_axis)
problem.SetParameterBlockConstant(rot_vec);
if (data.refine_scale)
problem.SetParameterLowerBound(&data.scale_factor, 0, 0.0);
else
problem.SetParameterBlockConstant(&data.scale_factor);
if (!data.refine_B)
problem.SetParameterBlockConstant(&data.B_factor);
if (data.refine_R) {
problem.SetParameterLowerBound(data.R, 0, 1e-5);
problem.SetParameterLowerBound(data.R, 1, 1e-5);
} else {
problem.SetParameterBlockConstant(data.R);
}
// ---- Solve ----------------------------------------------------------------
ceres::Solver::Options options;
options.linear_solver_type = ceres::DENSE_QR;
options.minimizer_progress_to_stdout = false;
options.logging_type = ceres::LoggingType::SILENT;
options.max_solver_time_in_seconds = data.max_time_s;
options.num_threads = 1;
ceres::Solver::Summary summary;
ceres::Solve(options, &problem, &summary);
data.final_cost = summary.final_cost;
data.solved = summary.IsSolutionUsable();
// ---- Write back refined geometry + lattice --------------------------------
if (data.refine_beam_center)
data.geom.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]);
if (data.refine_distance)
data.geom.DetectorDistance_mm(dist_mm);
if (data.refine_detector_angles)
data.geom.PoniRot1_rad(detector_rot[0]).PoniRot2_rad(detector_rot[1]);
if (data.refine_orientation || data.refine_unit_cell) {
double beta = data.latt.GetUnitCell().beta;
switch (data.crystal_system) {
case gemmi::CrystalSystem::Orthorhombic:
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, M_PI/2, M_PI/2);
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
break;
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);
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
latt_vec1[0] = (latt_vec1[0] + latt_vec1[1]) / 2.0;
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);
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
latt_vec1[0] = (latt_vec1[0] + latt_vec1[1] + latt_vec1[2]) / 3.0;
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);
LatticeToRodriguesAndLengths_Hex(data.latt, latt_vec0, latt_vec1);
break;
case gemmi::CrystalSystem::Monoclinic:
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, latt_vec2[0], M_PI/2);
LatticeToRodriguesLengthsBeta_Mono(data.latt, latt_vec0, latt_vec1, beta);
latt_vec2[0] = beta;
break;
default:
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1,
latt_vec2[0], latt_vec2[1], latt_vec2[2]);
default: {
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
const auto uc = data.latt.GetUnitCell();
latt_vec2[0] = uc.alpha * M_PI / 180.0;
latt_vec2[1] = uc.beta * M_PI / 180.0;
latt_vec2[2] = uc.gamma * M_PI / 180.0;
break;
}
}
}
// ---- 4. Build the problem ---------------------------------------------
ceres::Problem problem;
for (const auto &g : groups) {
for (const auto &obs : g.pixels) {
auto *cost = new ceres::AutoDiffCostFunction<
PixelResidual, 1, 2, 1, 2, 3, 3, 3, 3, 1, 1, 2>(
new PixelResidual(obs, g.Itrue, lambda, pixel_size,
g.h, g.k, g.l, data.crystal_system));
problem.AddResidualBlock(cost, new ceres::HuberLoss(3.0),
beam, &dist_mm, detector_rot, rot_vec,
latt_vec0, latt_vec1, latt_vec2,
&data.scale_factor, &data.B_factor, data.R);
}
}
data.residual_count = problem.NumResidualBlocks();
// ---- 5. Constrain / bound parameter blocks ----------------------------
if (!data.refine_orientation)
problem.SetParameterBlockConstant(latt_vec0);
if (!data.refine_unit_cell) {
problem.SetParameterBlockConstant(latt_vec1);
problem.SetParameterBlockConstant(latt_vec2);
} else {
for (int i = 0; i < 3; ++i) {
problem.SetParameterLowerBound(latt_vec1, i, 5.0);
problem.SetParameterUpperBound(latt_vec1, i, 1000.0);
}
if (data.crystal_system != gemmi::CrystalSystem::Monoclinic &&
data.crystal_system != gemmi::CrystalSystem::Triclinic)
problem.SetParameterBlockConstant(latt_vec2);
}
if (!data.refine_beam_center)
problem.SetParameterBlockConstant(beam);
if (!data.refine_distance) {
problem.SetParameterBlockConstant(&dist_mm);
} else {
problem.SetParameterLowerBound(&dist_mm, 0, dist_mm * 0.9);
problem.SetParameterUpperBound(&dist_mm, 0, dist_mm * 1.1);
}
if (!data.refine_detector_angles) {
problem.SetParameterBlockConstant(detector_rot);
} else {
const double rng = 3.0 / 180.0 * M_PI;
for (int i = 0; i < 2; ++i) {
problem.SetParameterLowerBound(detector_rot, i, detector_rot[i] - rng);
problem.SetParameterUpperBound(detector_rot, i, detector_rot[i] + rng);
}
}
if (!data.refine_rotation_axis)
problem.SetParameterBlockConstant(rot_vec);
if (data.refine_scale)
problem.SetParameterLowerBound(&data.scale_factor, 0, 0.0);
else
problem.SetParameterBlockConstant(&data.scale_factor);
if (!data.refine_B)
problem.SetParameterBlockConstant(&data.B_factor);
if (data.refine_R) {
problem.SetParameterLowerBound(data.R, 0, 1e-5);
problem.SetParameterLowerBound(data.R, 1, 1e-5);
} else {
problem.SetParameterBlockConstant(data.R);
}
// ---- 6. Solve ---------------------------------------------------------
ceres::Solver::Options options;
options.linear_solver_type = ceres::DENSE_QR;
options.minimizer_progress_to_stdout = false;
options.logging_type = ceres::LoggingType::SILENT;
options.max_solver_time_in_seconds = data.max_time_s;
options.num_threads = 1;
ceres::Solver::Summary summary;
ceres::Solve(options, &problem, &summary);
data.final_cost = summary.final_cost;
data.solved = summary.IsSolutionUsable();
// ---- 7. Write refined geometry + lattice back into data ---------------
if (data.refine_beam_center)
data.geom.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]);
if (data.refine_distance)
data.geom.DetectorDistance_mm(dist_mm);
if (data.refine_detector_angles)
data.geom.PoniRot1_rad(detector_rot[0]).PoniRot2_rad(detector_rot[1]);
if (data.refine_orientation || data.refine_unit_cell) {
switch (data.crystal_system) {
case gemmi::CrystalSystem::Orthorhombic:
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, M_PI/2, M_PI/2);
break;
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;
default:
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1,
latt_vec2[0], latt_vec2[1], latt_vec2[2]);
break;
}
}
} // predict<->refine iterations
// ---- Extract integrated reflections ---------------------------------------
// Profile-fitted intensity with the optimized profile P (shape only):
// I = sum_p [ P_p (Iobs_p - Ibkg_p) / v_p ] / sum_p [ P_p^2 / v_p ]
// var(I) = 1 / sum_p [ P_p^2 / v_p ]
// where P_p is the *shape* of the model per pixel (signal at G=Itrue=B=1).
// Two quantities are read back per reflection, using the *optimized* model:
//
// * Recorded (partial) intensity J by profile fitting against the normalized
// tangential profile P_t (sum over the shoebox ~ 1):
// J = sum_p [ P_t,p (Iobs_p - Ibkg_p) / v_p ] / sum_p [ P_t,p^2 / v_p ]
// var(J) = 1 / sum_p [ P_t,p^2 / v_p ]
// Because P_t is area-normalized, J estimates the integrated intensity
// actually recorded on this image (not the full-spot intensity).
//
// * Partiality: the profile-weighted mean of the *peak-normalized* radial
// factor exp(-eps_r^2/R0^2) in (0,1], i.e. how close to the Ewald sphere
// the reflection sits. Kept dimensionless for consistency with the rest of
// the pipeline (image_scale_corr = rlp / partiality, full I = J / partiality).
data.reflections.reserve(groups.size());
for (const auto &g : groups) {
double num = 0.0, den = 0.0, partiality_sum = 0.0, bkg_sum = 0.0;
double num = 0.0, den = 0.0, bkg_sum = 0.0;
double radial_sum = 0.0, radial_w = 0.0;
size_t n = 0;
for (const auto &obs : g.pixels) {
@@ -585,16 +633,19 @@ void PixelRefine::Run(const T *image,
if (!(data.R[0] > 0.0) || !(data.R[1] > 0.0))
continue;
// Normalized tangential profile (sum over shoebox ~ 1) -> fit weight.
const double P_t = obs.A_recip * std::exp(-eps_t_sq / (data.R[1] * data.R[1]))
/ (M_PI * data.R[1] * data.R[1]);
// Peak-normalized radial factor (the partiality), in (0,1].
const double P_radial = std::exp(-eps_r * eps_r / (data.R[0] * data.R[0]));
const double norm_tang = obs.A_recip / (M_PI * data.R[1] * data.R[1]);
const double P_tang = norm_tang * std::exp(-eps_t_sq / (data.R[1] * data.R[1]));
const double P = P_radial * P_tang; // model shape per pixel (unit scale)
const double v = SafeInv(obs.weight * obs.weight, 1.0); // pixel variance
const double signal = obs.Iobs - obs.Ibkg;
num += P * signal / v;
den += P * P / v;
partiality_sum += P_tang; // tangential profile sums to ~partiality
num += P_t * signal / v;
den += P_t * P_t / v;
radial_sum += P_radial * P_t; // weight partiality by the spot core
radial_w += P_t;
bkg_sum += obs.Ibkg;
++n;
}
@@ -609,11 +660,10 @@ void PixelRefine::Run(const T *image,
r.observed_x = NAN;
r.observed_y = NAN;
r.rlp = 1.0f;
r.partiality = static_cast<float>(partiality_sum);
r.partiality = (radial_w > 0.0) ? static_cast<float>(radial_sum / radial_w) : 1.0f;
if (den > 0.0 && n > 0) {
const double I = num / den;
r.I = static_cast<float>(I);
r.I = static_cast<float>(num / den);
r.sigma = static_cast<float>(std::sqrt(1.0 / den));
r.bkg = static_cast<float>(bkg_sum / static_cast<double>(n));
r.observed = true;