diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index 8f78e16c..225778f7 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -993,6 +993,122 @@ std::vector PixelRefine::PredictImage(const AzimuthalIntegrationProfile & return img; } +template +std::vector PixelRefine::ChiSquaredImage(const T *image, + const AzimuthalIntegrationProfile &profile, + BraggPrediction &prediction, + const PixelRefineData &data) const { + std::vector img(xpixel * ypixel, 0.0f); + + const double lambda = data.geom.GetWavelength_A(); + const double pixel_size = data.geom.GetPixelSize_mm(); + const auto azim_result = profile.GetResult(); + const auto azim_std = profile.GetStd(); + const auto &pixel_to_bin = mapping.GetPixelToBin(); + const auto &corrections = mapping.Corrections(); + const int total_bin_count = static_cast(azim_result.size()); + const double angle_rad = data.angle_deg * M_PI / 180.0; + const int radius = data.shoebox_radius; + const double bw = data.bandwidth; + + 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(); + }; + auto bandwidth_radial_sq = [&](double d) -> double { + if (bw <= 0.0 || d <= 0.0) + return 0.0; + const double bl = bw * lambda; + return bl * bl / (2.0 * d * d * d * d); + }; + + double beam[2], dist_mm, detector_rot[2], rot_vec[3]; + double latt_vec0[3], latt_vec1[3], latt_vec2[3]; + BuildParameterBlocks(data, beam, dist_mm, detector_rot, rot_vec, latt_vec0, latt_vec1, latt_vec2); + + 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 BraggPredictionSettings settings_prediction{ + .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(), + .max_hkl = 100, + .centering = data.centering, + .bandwidth_sigma = static_cast(data.bandwidth) + }; + const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction); + const auto &predicted = prediction.GetReflections(); + + for (int ri = 0; ri < nrefl; ++ri) { + const auto &refl = predicted[ri]; + const auto it = reference_data.find(hkl_key_generator(refl)); + if (it == reference_data.end()) + continue; + + const double Itrue = it->second; + const double R_bw_sq = bandwidth_radial_sq(refl.d); + + 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]; + + // Same gating as Run(): only pixels that actually enter the fit. + if (azim_bin >= total_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]; + const double Ibkg_sigma = azim_std[azim_bin]; + const double raw = static_cast(image[npixel]); + const double Iobs = raw * correction; + + double var = correction * std::max(Iobs, 0.0) + Ibkg_sigma * Ibkg_sigma; + if (!(var > 1.0)) + var = 1.0; + const double weight = 1.0 / std::sqrt(var); + + PixelObs obs{ + .x = static_cast(x), + .y = static_cast(y), + .Iobs = Iobs, + .Ibkg = Ibkg, + .weight = weight, + .A_recip = recip_area(x, y), + .angle_rad = angle_rad + }; + PixelResidual pr(obs, Itrue, lambda, pixel_size, + refl.h, refl.k, refl.l, R_bw_sq, data.crystal_system); + + double Ipred = 0.0; + if (pr.Model(beam, &dist_mm, detector_rot, rot_vec, + latt_vec0, latt_vec1, latt_vec2, + &data.scale_factor, &data.B_factor, data.R, Ipred)) { + // residual_i = (I_pred - I_obs) * weight (== Ceres residual); + // its square is this pixel's contribution to the cost. + const double rw = (Ipred - Iobs) * weight; + img[npixel] += static_cast(rw * rw); + } + } + } + } + + return img; +} + // Explicit instantiations for the supported (uncompressed) image pixel types. template void PixelRefine::Run(const int8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); template void PixelRefine::Run(const int16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); @@ -1000,3 +1116,10 @@ template void PixelRefine::Run(const int32_t *, const AzimuthalIntegrat template void PixelRefine::Run(const uint8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); template void PixelRefine::Run(const uint16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); template void PixelRefine::Run(const uint32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); + +template std::vector PixelRefine::ChiSquaredImage(const int8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::ChiSquaredImage(const int16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::ChiSquaredImage(const int32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::ChiSquaredImage(const uint8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::ChiSquaredImage(const uint16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; +template std::vector PixelRefine::ChiSquaredImage(const uint32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, const PixelRefineData &) const; diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index 9522fc91..73a63094 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -182,4 +182,16 @@ public: BraggPrediction &prediction, const PixelRefineData &data, bool include_background = true) const; + + // Render the per-pixel chi-square (cost density) that the optimizer actually + // minimizes: for every shoebox pixel that enters the fit it stores the squared + // weighted residual ((I_pred - I_obs)/sigma)^2 in *corrected* units - identical + // to the Ceres residual_i^2 - accumulating where shoeboxes overlap. Pixels that + // are not part of any shoebox stay 0; masked/saturated pixels (skipped by the + // fit) also stay 0. Summing the image gives ~2*final_cost. Diagnostic tool. + template + std::vector ChiSquaredImage(const T *image, + const AzimuthalIntegrationProfile &profile, + BraggPrediction &prediction, + const PixelRefineData &data) const; }; diff --git a/viewer/JFJochImageReadingWorker.cpp b/viewer/JFJochImageReadingWorker.cpp index 6582700c..8f90cb11 100644 --- a/viewer/JFJochImageReadingWorker.cpp +++ b/viewer/JFJochImageReadingWorker.cpp @@ -8,6 +8,7 @@ #include #include "JFJochImageReadingWorker.h" +#include "../reader/JFJochReaderImage.h" // JFJochReaderImage + GAP/ERROR/SATURATED sentinels #include "../image_analysis/LoadFCalcFromMtz.h" #include "../image_analysis/bragg_prediction/BraggPredictionFactory.h" #include "../image_analysis/geom_refinement/AssignSpotsToRings.h" @@ -69,6 +70,7 @@ JFJochImageReadingWorker::JFJochImageReadingWorker(const SpotFindingSettings &se indexing_settings(experiment.GetIndexingSettings()), azint_settings(experiment.GetAzimuthalIntegrationSettings()) { qRegisterMetaType("PixelRefineParams"); + qRegisterMetaType>("QVector"); spot_finding_settings = settings;; indexing = std::make_unique(indexing_settings); @@ -798,6 +800,74 @@ std::shared_ptr JFJochImageReadingWorker::WrapFloatImage_i(const st return si; } +void JFJochImageReadingWorker::SquaredResidualWithImage_i(std::vector &pred) const { + // PredictImage() returns raw detector units (same as the measured counts), so + // pred - measured is the per-pixel residual the model fails to explain. We plot + // |pred - measured|^2: sign-free, so it needs no diverging colour scale and just + // highlights where the model disagrees most. Masked / saturated pixels carry + // sentinels rather than counts, so no comparison is possible -> NaN (gap). + if (!current_image_ptr) + return; + const auto &img = current_image_ptr->Image(); + const size_t n = std::min(pred.size(), img.size()); + for (size_t i = 0; i < n; ++i) { + const int32_t v = img[i]; + if (v == GAP_PXL_VALUE || v == ERROR_PXL_VALUE || v == SATURATED_PXL_VALUE) { + pred[i] = NAN; + } else { + const float diff = pred[i] - static_cast(v); + pred[i] = diff * diff; + } + } +} + +void JFJochImageReadingWorker::MaskMeasuredSentinels_i(std::vector &img) const { + // The chi^2 image is 0 outside shoeboxes; show masked/saturated pixels as a gap + // (NaN) instead, so they read as "not comparable" rather than "zero cost". + if (!current_image_ptr) + return; + const auto &measured = current_image_ptr->Image(); + const size_t n = std::min(img.size(), measured.size()); + for (size_t i = 0; i < n; ++i) { + const int32_t v = measured[i]; + if (v == GAP_PXL_VALUE || v == ERROR_PXL_VALUE || v == SATURATED_PXL_VALUE) + img[i] = NAN; + } +} + +QVector JFJochImageReadingWorker::BuildShoeboxes_i(const PixelRefineData &data) const { + // One rectangle per fitted reflection: the shoebox the optimizer summed over, + // centred on the predicted position with half-size data.shoebox_radius. + QVector boxes; + boxes.reserve(static_cast(data.reflections.size())); + const int r = data.shoebox_radius; + const int side = 2 * r + 1; + for (const auto &refl : data.reflections) { + if (!std::isfinite(refl.predicted_x) || !std::isfinite(refl.predicted_y)) + continue; + const int cx = static_cast(std::lround(refl.predicted_x)); + const int cy = static_cast(std::lround(refl.predicted_y)); + boxes.push_back(QRect(cx - r, cy - r, side, side)); + } + return boxes; +} + +std::vector JFJochImageReadingWorker::BuildDisplayImage_i(const PixelRefineData &data, + int display_mode) const { + if (display_mode == PixelRefineParams::ChiSquared) { + // The cost density the optimizer actually minimizes (weighted residual^2). + const auto &img32 = current_image_ptr->Image(); + auto chi2 = pixel_refine_->ChiSquaredImage(img32.data(), *last_profile_, *pixel_pred_, data); + MaskMeasuredSentinels_i(chi2); + return chi2; + } + + auto pred = pixel_refine_->PredictImage(*last_profile_, *pixel_pred_, data, true); + if (display_mode == PixelRefineParams::SquaredDifference) + SquaredResidualWithImage_i(pred); + return pred; +} + void JFJochImageReadingWorker::LoadReference(QString path) { QMutexLocker ul(&m); try { @@ -840,8 +910,9 @@ void JFJochImageReadingWorker::PixelRefinePreview(PixelRefineParams params) { pixel_refine_->Run(img32.data(), *last_profile_, *pixel_pred_, d); 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)); + auto display = BuildDisplayImage_i(d, params.display_mode); + emit predictedImageReady(WrapFloatImage_i(display)); + emit predictedShoeboxes(BuildShoeboxes_i(d)); } catch (const std::exception &e) { emit pixelRefineStatus(QString("PixelRefine preview failed: %1").arg(e.what())); } @@ -884,8 +955,9 @@ void JFJochImageReadingWorker::PixelRefineRun(PixelRefineParams params) { emit pixelRefineParamsRefined(out); 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)); + auto display = BuildDisplayImage_i(d, params.display_mode); + emit predictedImageReady(WrapFloatImage_i(display)); + emit predictedShoeboxes(BuildShoeboxes_i(d)); // Show the refined predictions on the main image too. auto new_image = std::make_shared(*current_image_ptr); diff --git a/viewer/JFJochImageReadingWorker.h b/viewer/JFJochImageReadingWorker.h index 93d4240b..501b557d 100644 --- a/viewer/JFJochImageReadingWorker.h +++ b/viewer/JFJochImageReadingWorker.h @@ -71,6 +71,17 @@ private: void EnsurePixelRefine_i(); bool BuildPixelSeed_i(PixelRefineData &d, const PixelRefineParams &p, QString &reason) const; std::shared_ptr WrapFloatImage_i(const std::vector &img) const; + // Turn a predicted image into the squared residual |predicted - measured|^2 in + // place. Masked/saturated pixels become NaN (rendered as a gap: no comparison + // possible), not 0. + void SquaredResidualWithImage_i(std::vector &pred) const; + // Mark masked/saturated pixels of the current image as NaN (gap) in a float + // image, leaving the rest untouched (used for the chi^2 view). + void MaskMeasuredSentinels_i(std::vector &img) const; + // Build the per-reflection shoebox rectangles for the last refine/preview. + QVector BuildShoeboxes_i(const PixelRefineData &data) const; + // Build the float image to display for the given PixelRefineParams::DisplayMode. + std::vector BuildDisplayImage_i(const PixelRefineData &data, int display_mode) const; std::unique_ptr roi; @@ -134,6 +145,7 @@ signals: // PixelRefine (experimental) void predictedImageReady(std::shared_ptr image); + void predictedShoeboxes(QVector boxes); // per-reflection optimization windows void pixelRefineResidual(double cost, double cc, int64_t n_reflections); void pixelRefineParamsRefined(PixelRefineParams params); void pixelRefineStatus(QString message); diff --git a/viewer/JFJochViewerWindow.cpp b/viewer/JFJochViewerWindow.cpp index de383008..9fa6a302 100644 --- a/viewer/JFJochViewerWindow.cpp +++ b/viewer/JFJochViewerWindow.cpp @@ -28,6 +28,7 @@ #include "windows/JFJochPixelRefineWindow.h" #include "windows/JFJochMagnifierWindow.h" #include "image_viewer/JFJochImage.h" +#include "image_viewer/JFJochSimpleImage.h" #include JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString &file) : QMainWindow(parent) { @@ -351,6 +352,8 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString reading_worker, &JFJochImageReadingWorker::LoadReference); connect(reading_worker, &JFJochImageReadingWorker::predictedImageReady, pixelRefineWindow, &JFJochPixelRefineWindow::setPredictedImage); + connect(reading_worker, &JFJochImageReadingWorker::predictedShoeboxes, + pixelRefineWindow->imageView(), &JFJochSimpleImage::setShoeboxes); connect(reading_worker, &JFJochImageReadingWorker::pixelRefineResidual, pixelRefineWindow, &JFJochPixelRefineWindow::setResidual); connect(reading_worker, &JFJochImageReadingWorker::pixelRefineParamsRefined, diff --git a/viewer/image_viewer/JFJochSimpleImage.cpp b/viewer/image_viewer/JFJochSimpleImage.cpp index 24bdcad6..d3494a8f 100644 --- a/viewer/image_viewer/JFJochSimpleImage.cpp +++ b/viewer/image_viewer/JFJochSimpleImage.cpp @@ -40,6 +40,31 @@ void JFJochSimpleImage::setImage(std::shared_ptr img) { } } +void JFJochSimpleImage::setShoeboxes(QVector boxes) { + shoeboxes_ = std::move(boxes); + // Redraw overlays on the current image (no-op if no image yet). + updateOverlay(); +} + +void JFJochSimpleImage::addCustomOverlay() { + if (shoeboxes_.isEmpty() || !scene()) + return; + + // Cosmetic 1-px outline so the box edges stay thin at any zoom; only draw the + // ones currently in view (there can be hundreds of reflections). + const QRectF visibleRect = mapToScene(viewport()->geometry()).boundingRect(); + QPen pen(QColor(0, 220, 255), 0); // cyan, distinct from the prediction colours + pen.setCosmetic(true); + + for (const QRect &b : shoeboxes_) { + const QRectF r(b.x(), b.y(), b.width(), b.height()); + if (!visibleRect.intersects(r)) + continue; + auto *item = scene()->addRect(r, pen); + addOverlayItem(item); + } +} + void JFJochSimpleImage::mouseHover(QMouseEvent *event) { if (image_) { const QPointF scenePos = mapToScene(event->pos()); diff --git a/viewer/image_viewer/JFJochSimpleImage.h b/viewer/image_viewer/JFJochSimpleImage.h index 8d7251e4..40abac06 100644 --- a/viewer/image_viewer/JFJochSimpleImage.h +++ b/viewer/image_viewer/JFJochSimpleImage.h @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -18,14 +19,21 @@ class JFJochSimpleImage : public JFJochImage { Q_OBJECT std::shared_ptr image_; + + // Per-reflection shoebox rectangles (pixel coordinates) to overlay: the pixels + // PixelRefine actually summed over. Empty = nothing drawn. + QVector shoeboxes_; + // Prepare image template void loadImageInternal(const uint8_t *input); void loadImageInternal(); void mouseHover(QMouseEvent *event) override; + void addCustomOverlay() override; public: explicit JFJochSimpleImage(QWidget *parent = nullptr); public slots: void setImage(std::shared_ptr img); + void setShoeboxes(QVector boxes); }; diff --git a/viewer/windows/JFJochPixelRefineWindow.cpp b/viewer/windows/JFJochPixelRefineWindow.cpp index 230d4974..7f7a64ec 100644 --- a/viewer/windows/JFJochPixelRefineWindow.cpp +++ b/viewer/windows/JFJochPixelRefineWindow.cpp @@ -31,6 +31,15 @@ JFJochPixelRefineWindow::JFJochPixelRefineWindow(QWidget *parent) auto controlsLayout = new QVBoxLayout(controls); layout->addWidget(controls, 0); + // --- what the left image shows ------------------------------------------ + m_displayMode = new QComboBox(this); + m_displayMode->addItem(tr("Prediction")); + m_displayMode->addItem(tr("Squared difference |pred - image|²")); + m_displayMode->addItem(tr("χ² (weighted residual = LSQ cost)")); + auto displayForm = new QFormLayout(); + displayForm->addRow(tr("Display:"), m_displayMode); + controlsLayout->addLayout(displayForm); + auto paramBox = new QGroupBox(tr("Model parameters"), this); auto form = new QFormLayout(paramBox); @@ -97,6 +106,8 @@ JFJochPixelRefineWindow::JFJochPixelRefineWindow(QWidget *parent) for (auto *s : {m_R0, m_R1, m_bw, m_scale, m_B, m_beamx, m_beamy}) connect(s, &SliderPlusBox::valueChanged, this, [this](double) { onControlChanged(); }); + connect(m_displayMode, &QComboBox::currentIndexChanged, this, [this](int) { onControlChanged(); }); + connect(m_overrideBeam, &QCheckBox::toggled, this, [this](bool on) { m_beamx->setEnabled(on); m_beamy->setEnabled(on); @@ -141,6 +152,7 @@ PixelRefineParams JFJochPixelRefineWindow::currentParams() const { p.refine_scale = m_refScale->isChecked(); p.refine_B = m_refB->isChecked(); p.refine_R = m_refR->isChecked(); + p.display_mode = m_displayMode->currentIndex(); return p; } diff --git a/viewer/windows/JFJochPixelRefineWindow.h b/viewer/windows/JFJochPixelRefineWindow.h index 7940ce16..99c628c2 100644 --- a/viewer/windows/JFJochPixelRefineWindow.h +++ b/viewer/windows/JFJochPixelRefineWindow.h @@ -10,6 +10,7 @@ #include #include +#include #include #include #include @@ -33,6 +34,8 @@ class JFJochPixelRefineWindow : public JFJochHelperWindow { SliderPlusBox *m_beamx; SliderPlusBox *m_beamy; + QComboBox *m_displayMode; // Prediction vs. Difference (prediction - image) + QCheckBox *m_overrideBeam; QCheckBox *m_refOrientation; QCheckBox *m_refCell; diff --git a/viewer/windows/PixelRefineParams.h b/viewer/windows/PixelRefineParams.h index 4b266e62..64f6e27b 100644 --- a/viewer/windows/PixelRefineParams.h +++ b/viewer/windows/PixelRefineParams.h @@ -28,6 +28,14 @@ struct PixelRefineParams { bool refine_R = true; int max_iterations = 3; // <=0 means evaluate-only (preview / residual) + + // Display only (no effect on the fit): what the preview/refine image shows. + enum DisplayMode : int { + Prediction = 0, // forward-model image + SquaredDifference = 1, // |prediction - measured|^2 (raw, unweighted) + ChiSquared = 2 // ((prediction - measured)/sigma)^2 = the LSQ cost density + }; + int display_mode = Prediction; }; Q_DECLARE_METATYPE(PixelRefineParams)