1.0.0-rc.174 (#84)

* Rugnux: Performance improvements on GPU and CPU (more of the pre-scan and of scaling on the GPU, faster CPU spot finding and crystal refinement), with unchanged results.
* Rugnux: More robust processing - patches of persistently hot pixels are masked, an inconsistent merge triggers a retry at the measured beam centre, and builds targeting different CPU levels give the same results.
* Rugnux: Improved scaling and merging - reflections with an overloaded pixel are dropped, as in XDS, sparse rotation sweeps are scaled more reliably, and French-Wilson amplitudes use an anisotropic Wilson prior.
* Rugnux: Improved space-group determination - glide planes in groups without a centre of symmetry, screw axes from short or weak axial rows kept when a higher group is adopted, and more reliable decisions on twinned and pseudo-symmetric crystals.
* Rugnux: Improved small-molecule processing - spots that grow wider than the integration disk and split spots are integrated over their measured footprint, sparse lattices are integrated on every frame, and the `.hkl` file holds unmerged scaled reflections (SHELX HKLF 4).
* Rugnux: Reads Rigaku d*TREK SMV images (Saturn CCD), including detector 2theta and encoded pixel overflows; home-source (rotating-anode) datasets were added to the validation battery.
* jfjoch_viewer: Fixed processing failing at the end with "Wrong JPEG library version" on Linux; the merge window shows the space group with proper subscripts and a checklist of crystal pathologies.

Reviewed-on: #84
Co-authored-by: Filip Leonarski <filip.leonarski@psi.ch>

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -62,8 +62,10 @@ struct ImageSampler {
 
 }   // namespace
 
-BraggIntegrationEngineCPU::BraggIntegrationEngineCPU(const DiffractionExperiment &experiment)
+BraggIntegrationEngineCPU::BraggIntegrationEngineCPU(const DiffractionExperiment &experiment,
+                                                     const PixelMask &mask)
     : BraggIntegrationEngine(experiment),
+      static_mask(mask.GetPackedMask()),
       refl_mask(static_cast<int>(xpixel), static_cast<int>(ypixel), 0),
       owner(static_cast<int>(xpixel), static_cast<int>(ypixel), BRAGG_OWNER_NONE) {}
 
@@ -158,6 +160,14 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         int cx = 0, cy = 0, shell = -1;
         bool ok = false, strong = false, has_obs = false;
         bool full = false;      // every pixel of the signal disk was readable
+        bool overloaded = false;   // a pixel of the signal disk was saturated
+    };
+    // A saturated pixel, or one unreadable on this frame alone - the overload marker of a detector
+    // that writes its error value for a pixel it could not count (EIGER), which the preprocessor turns
+    // into a masked pixel like a gap's. The run's mask holds the gaps, so whatever is unreadable beyond
+    // it was lost to the flux.
+    auto static_masked = [&](size_t idx) {
+        return idx / 32 < static_mask.size() && ((static_mask[idx / 32] >> (idx % 32)) & 1U);
     };
     std::vector<Rough> rough(npredicted);
     double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
@@ -224,7 +234,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
             for (int x = x0; x <= x1; ++x) {
                 const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
                 const int32_t px = img[y * W + x];
-                if (d.signal < r1_sq) {
+                if (BraggStencilInSignal(st, x - r.predicted_x, y - r.predicted_y, d.signal, r1_sq)) {
                     // A pixel a nearer neighbour owns carries that neighbour's flux, so in Exclude
                     // mode it leaves the disk entirely - the sum, the pixel count the background is
                     // subtracted with, and the all-or-nothing validity gate alike. The box sum only
@@ -235,6 +245,9 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                         else if (exclude) continue;
                     }
                     ++n_inner;
+                    if (px == INT32_MAX
+                        || (px == INT32_MIN && !static_masked(static_cast<size_t>(y) * W + x)))
+                        out.overloaded = true;
                     if (!valid(px)) continue;
                     I_sum += px;
                     I_sum_x += static_cast<int64_t>(x) * px;
@@ -269,7 +282,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         // reflection is dropped whole - the one thing the stencil geometry does to the DATA rather
         // than to a measurement. Counted here because it is the only direct evidence of a radius that has
         // outgrown the pattern it is integrating (BraggIntegrationCounts).
-        const bool keep_partial = full || mode != IntegratorMode::BoxSum;
+        const bool keep_partial = full || mode != IntegratorMode::BoxSum || out.overloaded;
         if (keep_partial && n_bkg <= 5) ++counts.bkg_starved;
         // Would every predicted neighbour, the tails included, leave the ring starved that the
         // detector alone would not? That is the pattern's density, what the guard on a widened
@@ -411,7 +424,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
             if (overlap == OverlapMode::Reject && rh.n_own < overlap_min_peak * rh.n_disk) continue;
             results[i] = {static_cast<float>(rh.I), static_cast<float>(rh.sigma), static_cast<float>(rh.bkg),
                           static_cast<float>(rh.obs_x), static_cast<float>(rh.obs_y),
-                          static_cast<float>(rh.var_bkg), true, rh.has_obs};
+                          static_cast<float>(rh.var_bkg), true, rh.has_obs, rh.overloaded};
         }
         return Finalize(predicted, npredicted, results, image_number);
     }
@@ -515,15 +528,26 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
     for (size_t i = 0; i < npredicted; ++i) {
         const auto &rh = rough[i];
         if (!rh.ok) continue;
+        // An overloaded reflection is not measured by any fit - its brightest pixels are what is
+        // missing. It is passed on with its box sum and the flag, for the merge to drop its rocking
+        // event (Reflection::overloaded); the missing-peak cut below would otherwise drop just this
+        // frame and leave the event to be extrapolated from its flanks.
+        if (rh.overloaded) {
+            results[i] = {static_cast<float>(rh.I), static_cast<float>(rh.sigma), static_cast<float>(rh.bkg),
+                          static_cast<float>(rh.obs_x), static_cast<float>(rh.obs_y),
+                          static_cast<float>(rh.var_bkg), true, false, true};
+            continue;
+        }
         const int sh = rh.shell < 0 ? 0 : rh.shell;
 
         int Rf = R;
         const std::vector<double> *Pvec = &shell_P[sh];
+        const BraggStencil st = MakeBraggStencil(predicted[i].predicted_x, predicted[i].predicted_y, stencil);
         if (!empirical) {
             const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y;
             const double Rpx = std::hypot(rx, ry);
             const double tan2t = Rpx / F_px;
-            const double s2t = shell_sigma2[sh].tan;
+            double s2t = shell_sigma2[sh].tan;
             double s2r = s2t, ux = 1.0, uy = 0.0;
             bool elong = false;
             if (use_ellipse) {
@@ -539,10 +563,18 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                     elong = true;
                 }
             }
+            // A spot that outgrows the r1 disk has widths the disk cannot measure - the moments above
+            // saturate near r1^2/4 - so the measured footprint sets them instead (SpotFootprint.h).
+            if (st.fp_s2t > 0.0f) {
+                ux = st.ux; uy = st.uy;
+                s2t = std::max<double>(s2t, st.fp_s2t);
+                s2r = std::max<double>(s2r, st.fp_s2r);
+                elong = true;
+            }
             // Build the Gaussian per reflection, centred on the sub-pixel predicted position and (when
-            // needed) radially elongated, on a grid grown to hold the streak.
+            // needed) elongated, on a grid grown to hold the wider of its two axes.
             const double fx = predicted[i].predicted_x - rh.cx, fy = predicted[i].predicted_y - rh.cy;
-            Rf = elong ? std::min(3 * R, static_cast<int>(std::ceil(r2 + 2.0 * std::sqrt(s2r)))) : R;
+            Rf = elong ? std::min(3 * R, static_cast<int>(std::ceil(r2 + 2.0 * std::sqrt(std::max(s2r, s2t))))) : R;
             // The grid's image rows, asked for now (as in pass A): building the profile below takes
             // long enough for them to arrive before the fit reads them.
             const int gx0 = std::max(0, rh.cx - Rf), gx1 = std::min(W - 1, rh.cx + Rf);
@@ -582,7 +614,14 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
                 if (Pp <= 0.0) continue;
                 const int x = rh.cx + dx, y = rh.cy + dy;
                 if (x < 0 || y < 0 || x >= W || y >= H) continue;
-                const bool in_disk = dx * dx + dy * dy < r1_sq;
+                // The region the summation seed was taken over (pass A), so the guard below
+                // compares the two over the same pixels.
+                const bool in_disk = st.fp_s2t > 0.0f
+                        ? BraggStencilInSignal(st, x - predicted[i].predicted_x, y - predicted[i].predicted_y,
+                                               (x - predicted[i].predicted_x) * (x - predicted[i].predicted_x)
+                                               + (y - predicted[i].predicted_y) * (y - predicted[i].predicted_y),
+                                               r1_sq)
+                        : dx * dx + dy * dy < r1_sq;
                 p_grid += Pp;
                 p_peak = std::max(p_peak, Pp);
                 if (in_disk) m_all += Pp;
@@ -652,7 +691,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
         results[i] = {static_cast<float>(I), static_cast<float>(sigma),
                       static_cast<float>(rh.bkg),
                       static_cast<float>(rh.obs_x), static_cast<float>(rh.obs_y),
-                      static_cast<float>(var_bkg), true, rh.has_obs};
+                      static_cast<float>(var_bkg), true, rh.has_obs, rh.overloaded};
     }
 
     return Finalize(predicted, npredicted, results, image_number);