diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 5dc02d55..c297de04 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -475,8 +475,13 @@ bool IndexAndRefine::PixelRefineIntegrate(DataMessage &msg, i_outcome.latt = prd.latt; i_outcome.image_scale_g = static_cast(prd.scale_factor); i_outcome.image_scale_b_factor_Ang2 = static_cast(prd.B_factor); + if (std::isfinite(prd.cc)) { + i_outcome.image_scale_cc = static_cast(prd.cc); + i_outcome.image_scale_cc_n = prd.cc_n; + } msg.image_scale_factor = static_cast(prd.scale_factor); + msg.image_scale_cc = i_outcome.image_scale_cc; if (prd.B_factor != 0.0) msg.image_scale_b_factor = static_cast(prd.B_factor); diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index 3b6ccdfe..b2d2578f 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -766,19 +766,20 @@ void PixelRefine::Run(const T *image, } // predict<->refine iterations // ---- Extract integrated reflections --------------------------------------- - // Two quantities are read back per reflection, using the *optimized* model: + // Profile fitting gives the recorded amplitude (against the normalized + // tangential profile P_t): + // J = sum_p[ P_t,p (Iobs_p - Ibkg_p)/v_p ] / sum_p[ P_t,p^2 / v_p ] + // ~ G * Itrue * B_term * partiality (recorded intensity) + // var(J) = 1 / sum_p[ P_t,p^2 / v_p ] // - // * 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). + // Output split (Merge multiplies r.I * image_scale_corr and weights by + // 1/(sigma*image_scale_corr)^2 - see Merge.cpp): + // r.I = J / (B_term * partiality) = G * Itrue (B/partiality corrected) + // r.sigma = sqrt(var(J)) / (B_term * partiality) + // r.partiality = profile-weighted peak radial factor in (0,1] (Merge filter only) + // r.image_scale_corr = 1/G (per-image scale ONLY) + // so r.I * image_scale_corr = Itrue. B and partiality live on the intensity, + // G lives on image_scale_corr - one clean meaning per field. data.reflections.reserve(groups.size()); for (const auto &g : groups) { double num = 0.0, den = 0.0, bkg_sum = 0.0; @@ -826,16 +827,22 @@ void PixelRefine::Run(const T *image, r.partiality = (radial_w > 0.0) ? static_cast(radial_sum / radial_w) : 1.0f; if (den > 0.0 && n > 0) { - r.I = static_cast(num / den); - r.sigma = static_cast(std::sqrt(1.0 / den)); + const double I_amp = num / den; // ~ G*Itrue*B_term*partiality + const double sigma_amp = std::sqrt(1.0 / den); + const double B_term = std::exp(-data.B_factor / (4.0 * g.d * g.d)); + const double corr = static_cast(r.partiality) * B_term; // B & partiality r.bkg = static_cast(bkg_sum / static_cast(n)); r.observed = true; - // Put I onto a common (merge-ready) scale: I * image_scale_corr is the - // full-spot, scale/DW-corrected intensity. Fold in the refined G, the - // per-reflection B_term and the partiality (rlp = 1 here). - const double B_term = std::exp(-data.B_factor / (4.0 * g.d * g.d)); - const double denom = static_cast(r.partiality) * data.scale_factor * B_term; - r.image_scale_corr = (denom > 0.0) ? static_cast(r.rlp / denom) : NAN; + + if (corr > 0.0 && data.scale_factor > 0.0) { + r.I = static_cast(I_amp / corr); // = G*Itrue + r.sigma = static_cast(sigma_amp / corr); + r.image_scale_corr = static_cast(1.0 / data.scale_factor); // G only + } else { + r.I = static_cast(I_amp); + r.sigma = static_cast(sigma_amp); + r.image_scale_corr = NAN; + } } else { r.I = 0.0f; r.sigma = NAN; @@ -844,6 +851,36 @@ void PixelRefine::Run(const T *image, } data.reflections.push_back(r); } + + // ---- Per-image CC vs reference (the half/ref correlation diagnostic) ------- + // Pearson CC of the scaled estimate (r.I * image_scale_corr = Itrue_est) + // against the reference intensities, over the matched reflections. + { + double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; + size_t cn = 0; + for (const auto &r : data.reflections) { + if (!r.observed || !std::isfinite(r.I) || !std::isfinite(r.image_scale_corr)) + continue; + const auto it = reference_data.find(hkl_key_generator(r)); + if (it == reference_data.end()) + continue; + const double x = static_cast(r.I) * static_cast(r.image_scale_corr); + const double y = it->second; + if (!std::isfinite(x) || !std::isfinite(y)) + continue; + sx += x; sy += y; sxx += x * x; syy += y * y; sxy += x * y; ++cn; + } + data.cc = NAN; + data.cc_n = static_cast(cn); + if (cn >= 3) { + const double nd = static_cast(cn); + const double cov = sxy - sx * sy / nd; + const double vx = sxx - sx * sx / nd; + const double vy = syy - sy * sy / nd; + if (vx > 0.0 && vy > 0.0) + data.cc = cov / std::sqrt(vx * vy); + } + } } std::vector PixelRefine::PredictImage(const AzimuthalIntegrationProfile &profile, diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index ae805f4e..9522fc91 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -139,6 +139,8 @@ struct PixelRefineData { bool solved = false; double final_cost = NAN; size_t residual_count = 0; + double cc = NAN; // per-image CC of scaled intensities vs reference + int64_t cc_n = 0; // number of reflections in the CC }; class PixelRefine { diff --git a/viewer/JFJochImageReadingWorker.cpp b/viewer/JFJochImageReadingWorker.cpp index 1ad60bd0..08eb7ad7 100644 --- a/viewer/JFJochImageReadingWorker.cpp +++ b/viewer/JFJochImageReadingWorker.cpp @@ -816,7 +816,7 @@ void JFJochImageReadingWorker::PixelRefinePreview(PixelRefineParams params) { try { const auto &img32 = current_image_ptr->Image(); pixel_refine_->Run(img32.data(), *last_profile_, *pixel_pred_, d); - emit pixelRefineResidual(d.final_cost, static_cast(d.reflections.size())); + emit pixelRefineResidual(d.final_cost, d.cc, static_cast(d.reflections.size())); auto pred = pixel_refine_->PredictImage(*last_profile_, *pixel_pred_, d, true); emit predictedImageReady(WrapFloatImage_i(pred)); @@ -859,7 +859,7 @@ void JFJochImageReadingWorker::PixelRefineRun(PixelRefineParams params) { out.beam_x = d.geom.GetBeamX_pxl(); out.beam_y = d.geom.GetBeamY_pxl(); emit pixelRefineParamsRefined(out); - emit pixelRefineResidual(d.final_cost, static_cast(d.reflections.size())); + emit pixelRefineResidual(d.final_cost, d.cc, static_cast(d.reflections.size())); auto pred = pixel_refine_->PredictImage(*last_profile_, *pixel_pred_, d, true); emit predictedImageReady(WrapFloatImage_i(pred)); diff --git a/viewer/JFJochImageReadingWorker.h b/viewer/JFJochImageReadingWorker.h index 6ab9e429..a5b219df 100644 --- a/viewer/JFJochImageReadingWorker.h +++ b/viewer/JFJochImageReadingWorker.h @@ -134,7 +134,7 @@ signals: // PixelRefine (experimental) void predictedImageReady(std::shared_ptr image); - void pixelRefineResidual(double cost, int64_t n_reflections); + void pixelRefineResidual(double cost, double cc, int64_t n_reflections); void pixelRefineParamsRefined(PixelRefineParams params); void pixelRefineStatus(QString message); diff --git a/viewer/windows/JFJochPixelRefineWindow.cpp b/viewer/windows/JFJochPixelRefineWindow.cpp index 364cf5df..6cce6a99 100644 --- a/viewer/windows/JFJochPixelRefineWindow.cpp +++ b/viewer/windows/JFJochPixelRefineWindow.cpp @@ -161,9 +161,12 @@ void JFJochPixelRefineWindow::setPredictedImage(std::shared_ptrsetImage(std::move(image)); } -void JFJochPixelRefineWindow::setResidual(double cost, int64_t n_reflections) { - m_residual->setText(tr("Residual: %1 (%2 reflections)") +void JFJochPixelRefineWindow::setResidual(double cost, double cc, int64_t n_reflections) { + const QString cc_str = std::isfinite(cc) ? QString::number(cc * 100.0, 'f', 1) + "%" + : QStringLiteral("—"); + m_residual->setText(tr("Residual: %1 CC: %2 (%3 reflections)") .arg(cost, 0, 'g', 6) + .arg(cc_str) .arg(n_reflections)); } diff --git a/viewer/windows/JFJochPixelRefineWindow.h b/viewer/windows/JFJochPixelRefineWindow.h index 7e7f167c..b508497e 100644 --- a/viewer/windows/JFJochPixelRefineWindow.h +++ b/viewer/windows/JFJochPixelRefineWindow.h @@ -67,7 +67,7 @@ signals: public slots: void setPredictedImage(std::shared_ptr image); - void setResidual(double cost, int64_t n_reflections); + void setResidual(double cost, double cc, int64_t n_reflections); void setRefinedParams(PixelRefineParams params); void setStatus(QString message); };