From 6f6098d5083bfd2d08c413f5c3068066b53e8212 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Mon, 8 Jun 2026 13:44:39 +0200 Subject: [PATCH] PixelRefine: Work in progress (Claude) --- .../pixel_refinement/PixelRefine.cpp | 536 ++++++++++-------- image_analysis/pixel_refinement/PixelRefine.h | 67 +++ 2 files changed, 360 insertions(+), 243 deletions(-) diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index e734516c..63971bf2 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -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 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(refl.predicted_y - radius, 0); - const int max_y = std::min(refl.predicted_y + radius, ypixel - 1); - const int min_x = std::max(refl.predicted_x - radius, 0); - const int max_x = std::min(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::max()) - continue; - if (std::is_signed_v && (image[npixel] == std::numeric_limits::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(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(x), - .y = static_cast(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(refl.predicted_y - radius, 0); + const int max_y = std::min(refl.predicted_y + radius, ypixel - 1); + const int min_x = std::max(refl.predicted_x - radius, 0); + const int max_x = std::min(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::max()) + continue; + if (std::is_signed_v && (image[npixel] == std::numeric_limits::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(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(x), + .y = static_cast(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(partiality_sum); + r.partiality = (radial_w > 0.0) ? static_cast(radial_sum / radial_w) : 1.0f; if (den > 0.0 && n > 0) { - const double I = num / den; - r.I = static_cast(I); + r.I = static_cast(num / den); r.sigma = static_cast(std::sqrt(1.0 / den)); r.bkg = static_cast(bkg_sum / static_cast(n)); r.observed = true; diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index 5848cf15..44ae21c8 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -9,6 +9,72 @@ #include "../common/AzimuthalIntegrationProfile.h" #include "../scale_merge/HKLKey.h" +// ============================================================================= +// PixelRefine — one optimization to rule geometry, integration and scaling +// ============================================================================= +// +// Intent +// ------ +// Classical crystallographic data processing is a one-way pipeline: +// +// spot finding -> indexing -> geometry refinement -> integration -> scaling -> merging +// +// Each stage consumes the previous stage's output and never talks back. The +// integrator trusts the refined geometry; the scaler trusts the integrated +// intensities; nothing downstream is ever allowed to correct an upstream +// parameter. Post-refinement and profile fitting were the field's partial +// answers to this: post-refinement lets merged intensities nudge per-image +// orientation/cell/mosaicity, and profile fitting lets a learned spot shape +// improve weak-reflection intensities. But both are narrow back-channels bolted +// onto a feed-forward pipeline — there is no end-to-end gradient that flows from +// the raw detector pixels all the way back to every parameter at once. +// +// PixelRefine is an experiment in doing the whole thing as a *single* least +// squares problem. We write down, for every pixel in a reflection's shoebox, the +// expected counts as an explicit forward model +// +// I_pred(pixel) = G * I_true * B_term * P_radial * P_tangential + I_bkg +// +// and let Ceres autodiff back-propagate the per-pixel residuals into ALL of: +// * detector geometry (beam centre, distance, tilt) +// * crystal orientation + unit cell +// * overall scale G and Debye-Waller B +// * the reciprocal-space spot widths R = (radial, tangential) +// simultaneously. Geometry refinement, profile-fitted integration and scaling +// then stop being separate stages: they are different parameters of one model, +// coupled through the same pixels, with full backpropagation between them. Once +// the model is differentiable end-to-end, things that used to need bespoke code +// — mosaicity refinement, profile fitting, partiality — fall out "for free" as +// extra parameters of the same forward model. +// +// How I_true enters +// ----------------- +// I_true is NOT refined here. It is a *fixed hypothesis* for the duration of a +// pass: the current best merged estimate of each reflection's full intensity. +// The intended outer loop is iterative, like EM / self-consistent field: +// +// repeat over the whole dataset: +// run PixelRefine on every image with the current I_true reference +// re-merge the resulting intensities -> new I_true +// until the reference stops changing +// +// So a single PixelRefine call answers "given this intensity hypothesis, what +// geometry/scale/profile best explains these pixels, and what intensities do I +// read back out?", and the dataset-level loop refines the hypothesis itself. +// +// Inner predict<->refine loop +// --------------------------- +// Within one image we also iterate (max_iterations): Bragg prediction places the +// shoeboxes, we refine, the refined geometry/cell feed the next prediction, etc. +// Initially the model geometry equals the experiment's DiffractionGeometry, but +// as refinement proceeds it diverges, so later predictions must use the *refined* +// data.geom rather than the static experiment geometry. +// +// Status: experimental prototype. The forward model (esp. the still-image +// Lorentz/partiality normalization) is deliberately simple and expected to +// evolve. See PixelRefine.cpp for the physics conventions and known caveats. +// ============================================================================= + struct PixelRefineData { // --- model state (input as initial guess, output as refined result) --- DiffractionGeometry geom; @@ -37,6 +103,7 @@ struct PixelRefineData { double max_time_s = 5.0; int shoebox_radius = 3; // half-size of the per-reflection pixel box + int max_iterations = 3; // inner predict<->refine cycles (re-predict with refined geom/latt) // --- output --- std::vector reflections; // profile-fitted integration result