diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index e0094f00..b1fcd14d 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -9,7 +9,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * Bragg integration: the profile fit's `background_variance` now takes the fitted intensity itself out of the fit variance instead of `max(0, I)`, so a reflection that fluctuated below zero no longer reports a background variance two to three times too small and is no longer weighted up for it. * Scaling: the rotation merge weights each combined full by its variance rebuilt at the reflection's mean intensity rather than by the full's own sigma, as the stills merge already did. * Scaling: a rotation frame whose fitted scale collapses - it recorded no diffraction from the indexed lattice - is now **dropped from the merge** instead of being merged unscaled, which had asserted a scale of 1 for a frame demonstrably nowhere near it. -* Bragg integration: a reflection whose signal disk is cut by a mask, an untrusted region, a detector gap or an overload is now profile-fitted over the pixels that remain instead of being discarded, as long as at least `--overlap-minpk` of its expected profile is readable (XDS's MINPK); `--integrator boxsum` still discards it. +* Bragg integration: a reflection whose signal disk is cut by a mask, an untrusted region, a detector gap or an overload is now profile-fitted over the pixels that remain instead of being discarded, as long as at least `--overlap-minpk` of its expected profile is readable (XDS's MINPK) **and the unreadable part does not take the profile's peak**; `--integrator boxsum` still discards it. * rugnux: Rotation data are integrated on **every frame whose spots the sweep's lattice explains**, instead of only on frames that would also index on their own; the reported indexing rate still counts the latter. * Scaling: a rotation frame too sparse to fit a rocking width of its own now takes the run's median instead of a fixed default. * rugnux: The error-model **a** and **b** are reported in XDS's convention, and `_reflns.jfjoch_diffrn_ISa` now carries the whole-range `1/sqrt(a*b)` that XDS's ISa denotes; the strong-reflection asymptote moves to `_reflns.jfjoch_diffrn_ISa_asymptotic`. **A file written by an earlier version carries the asymptote under the old name.** diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 90e55694..ee854bd0 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -640,6 +640,8 @@ where $c$ is the pixel value and the de-biased variance $v$ (background plus mod **Pixels the fit cannot use (MINPK).** A profile fit is the amplitude of a *normalised* profile, so a pixel left out of the sum renormalises the estimator by construction: it costs information — $\sum P^2/v$ shrinks and $\sigma$ grows — but biases nothing. That is what keeps a reflection whose signal disk is cut by a mask, an untrusted region, a detector gap or an overload: those pixels are simply not read, and the fit is taken over the rest, exactly as the shared pixels of a crowded reflection are (`--overlap exclude`). The reflection is kept only while enough of the expected profile survives — at least `--overlap-minpk` of the profile mass that falls on the detector at all, default 0.75, which is XDS's `MINPK` and dials' `valid_foreground_threshold`. The complete reflections alone teach the profile, its resolution shells and their widths. `--integrator boxsum` has no profile to renormalise with and keeps the all-or-nothing rule of §9.2. +"Biases nothing" holds while the pixels go missing for reasons that have nothing to do with the reflection, which is true of a gap, a mask or the edge of the sensor — where the loss is set by the detector, and the same hole recurs on every frame of the rocking curve because the spot does not move off it. It is not true of a pixel invalidated *by the flux it saw*: that pixel goes missing **because** the reflection was bright, and it is the peak. The fit then has only the wings to set the amplitude from and reads low — measured at $-50\%$ against the symmetry mates, on reflections carrying several times the mean intensity of their shell, which are the largest terms of $R_\mathrm{meas}$. MINPK cannot separate the two cases, because it cuts on profile *mass* and the peak of a broad spot is a few percent of the mass. So a second condition applies alongside it: **no unreadable pixel may carry more than 0.9 of the profile's own peak value**. As a fraction of the peak rather than a radius in pixels, that scales with the spot — for a Gaussian it is a cut at $\sqrt{-2\ln f}\,\sigma = 0.46\sigma$, the peak pixel alone where $\sigma$ is 0.8 px and the crest of the ridge where the profile is a bandwidth streak — and it needs nothing the fit does not already compute. It costs 0.05–0.07 % of the recovered observations. + The integrator is selected by `--integrator boxsum|gaussian|empirical` (default `gaussian`). ### 9.4 Lorentz–polarization factor handling diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 9b05f6ad..4226ac93 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -84,6 +84,22 @@ constexpr float MAX_STENCIL_GROW_OVER_R3 = 2.0f; // summation (box-sum) intensity when the profile result disagrees with the summation seed by more than // this many box-sum sigmas (a real fit agrees within counting noise, so the margin is generous). constexpr double PROFILE_SUMMATION_MAX_NSIGMA = 10.0; + +// MINPK keeps a reflection while enough of its expected profile is readable. It says nothing about +// WHERE the unreadable part is, and the two are not the same question. A pixel lost to a gap, a mask +// or the edge of the sensor is lost for reasons that have nothing to do with this reflection, and the +// fit renormalises over what is left with no bias. A pixel lost because the flux it saw put it over +// the detector's range is lost BECAUSE the reflection was bright, and it is the peak: the fit then has +// only the wings to set the amplitude from, and reads low - measured at -50% on the strongest +// low-resolution reflections, which are also the largest terms of R_meas. So no unreadable pixel may +// carry more than this fraction of the profile's own peak value. +// +// A fraction of the peak rather than a radius in pixels, because the peak is as wide as the spot: for +// a Gaussian the cut sits at sqrt(-2 ln f) sigma, i.e. 0.46 sigma here, which is the peak pixel alone +// where sigma is 0.8 px and the crest of the ridge where it is 2.4 px or a bandwidth streak. It also +// needs nothing the fit does not already have, so it costs one max-reduction in the loop that +// measures the readable fraction, and it applies to the learned empirical profile unchanged. +constexpr double MINPK_MAX_MISSING_PEAK = 0.9; } // namespace bragg_engine // One reflection's extracted intensity, produced by the derived engine and turned into a diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 18dc85ae..a6d243d6 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -479,6 +479,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, // summation seed the runaway guard compares against actually summed; with nothing missing // and nothing excluded it is 1 and the guard is untouched. --- double p_grid = 0.0, p_valid = 0.0, p_own = 0.0, m_all = 0.0, m_read = 0.0; + double p_peak = 0.0, p_lost_peak = 0.0; for (int dy = -Rf; dy <= Rf; ++dy) for (int dx = -Rf; dx <= Rf; ++dx) { const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)]; @@ -487,14 +488,23 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (x < 0 || y < 0 || x >= W || y >= H) continue; const bool in_disk = dx * dx + dy * dy < r1_sq; p_grid += Pp; + p_peak = std::max(p_peak, Pp); if (in_disk) m_all += Pp; - if (!valid(img[y * W + x])) continue; + if (!valid(img[y * W + x])) { + p_lost_peak = std::max(p_lost_peak, Pp); + continue; + } p_valid += Pp; const bool own = clean(x, y, i); if (own) p_own += Pp; if (in_disk && (own || !exclude)) m_read += Pp; } if (p_valid < overlap_min_peak * p_grid) continue; + // A hole in the profile's PEAK is a different defect from a hole in its wings, and the mass + // fraction above cannot tell them apart - the peak of a broad spot is a few percent of the + // mass, so MINPK passes a reflection that has lost the one part of the profile its amplitude + // is determined by. See the header. + if (p_lost_peak > MINPK_MAX_MISSING_PEAK * p_peak) continue; if (overlap == OverlapMode::Reject && p_own < overlap_min_peak) continue; const double B = std::max(rh.bkg, PIXEL_VARIANCE_FLOOR); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 178475c9..59d78686 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -34,6 +34,7 @@ struct BraggGpuParams { float claim_sq, inv_claim; // how far a reflection claims pixels in the owner map float minpk; // least readable/clean profile fraction that is kept (XDS MINPK) int partial_ok; // keep a signal disk with unreadable pixels: anything but a box sum + float peak_frac; // most of the profile's peak value an unreadable pixel may carry }; __device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; } @@ -467,6 +468,9 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x extern __shared__ float Pbuf[]; __shared__ float s_gs, s_num, s_den, s_I, s_wsum; __shared__ float s_pgrid, s_pvalid, s_pown, s_mall, s_mread; + // Max reductions. The profile is non-negative, and for non-negative floats the IEEE bit pattern + // orders exactly as the value does, so an integer atomicMax on that pattern is an exact float max. + __shared__ int s_ppeak_i, s_plost_i; __shared__ int s_Rf, s_Gf; if (!ok_a[i]) { if (threadIdx.x == 0) ok_o[i] = 0; return; } @@ -526,9 +530,11 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x // CPU engine. if (threadIdx.x == 0) { s_pgrid = 0.0f; s_pvalid = 0.0f; s_pown = 0.0f; s_mall = 0.0f; s_mread = 0.0f; + s_ppeak_i = 0; s_plost_i = 0; } __syncthreads(); float l_pgrid = 0.0f, l_pvalid = 0.0f, l_pown = 0.0f, l_mall = 0.0f, l_mread = 0.0f; + float l_ppeak = 0.0f, l_plost = 0.0f; for (int k = threadIdx.x; k < GfGf; k += blockDim.x) { const float Pp = Pbuf[k]; if (Pp <= 0.0f) continue; @@ -537,8 +543,12 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; const bool in_disk = (float) (dx * dx + dy * dy) < p.r1_sq; l_pgrid += Pp; + l_ppeak = fmaxf(l_ppeak, Pp); if (in_disk) l_mall += Pp; - if (!valid(img[y * p.W + x])) continue; + if (!valid(img[y * p.W + x])) { + l_plost = fmaxf(l_plost, Pp); + continue; + } l_pvalid += Pp; const bool own = !p.overlap || BraggOwnedBy(owner[y * p.W + x], i); // Zeroing the profile here is how Exclude drops the pixel: the fit skips any cell with @@ -550,11 +560,17 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x } atomicAdd(&s_pgrid, l_pgrid); atomicAdd(&s_pvalid, l_pvalid); atomicAdd(&s_pown, l_pown); atomicAdd(&s_mall, l_mall); atomicAdd(&s_mread, l_mread); + atomicMax(&s_ppeak_i, __float_as_int(l_ppeak)); atomicMax(&s_plost_i, __float_as_int(l_plost)); __syncthreads(); if (s_pvalid < p.minpk * s_pgrid) { if (threadIdx.x == 0) ok_o[i] = 0; return; } + // A hole in the profile's PEAK is a different defect from a hole in its wings. See the CPU engine. + if (__int_as_float(s_plost_i) > p.peak_frac * __int_as_float(s_ppeak_i)) { + if (threadIdx.x == 0) ok_o[i] = 0; + return; + } if (p.overlap == 1 && s_pown < p.minpk) { if (threadIdx.x == 0) ok_o[i] = 0; return; @@ -745,6 +761,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu .exclude_px = (overlap == OverlapMode::Exclude && mode != IntegratorMode::BoxSum) ? 1 : 0, .claim_sq = claim * claim, .inv_claim = inv_claim, .minpk = overlap_min_peak, .partial_ok = mode != IntegratorMode::BoxSum ? 1 : 0, + .peak_frac = static_cast(MINPK_MAX_MISSING_PEAK), }; // Whether the radial correction runs for THIS image. n_rad only says the buffers exist - under diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index e869db18..b0a0f899 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -41,7 +41,12 @@ Reflection MakeReflection(float x, float y, float d, int hkl) { // companion_dx > 0 puts a second spot that many pixels beside every grid spot, so their r1 signal // disks share pixels while the background rings still see clean sky - which is what a dense pattern // actually looks like (crowded along one reciprocal axis, sparse across it). -Scene BuildScene(size_t width, size_t height, int spacing = 60, float companion_dx = 0.0f) { +// clip_spots punches unreadable pixels into the spots themselves rather than into empty sky: the +// centre of every 5th, a mid-profile pixel of every 7th and a disk-edge pixel of every 11th. That is +// the MINPK rescue's own case - a reflection kept and fitted over the pixels it has - and with it the +// peak-loss rule, which has to fire on the same reflections in both engines. +Scene BuildScene(size_t width, size_t height, int spacing = 60, float companion_dx = 0.0f, + bool clip_spots = false) { Scene s; s.width = width; s.height = height; @@ -91,6 +96,19 @@ Scene BuildScene(size_t width, size_t height, int spacing = 60, float companion_ const size_t idx = (static_cast(k) * 2654435761u) % s.image.size(); s.image[idx] = (k % 2) ? INT32_MIN : INT32_MAX; } + + if (clip_spots) + for (size_t n = 0; n < s.predicted.size(); ++n) { + int dx = 0, dy = 0; + if (n % 5 == 0) { dx = 0; dy = 0; } // the peak itself: the rule must reject + else if (n % 7 == 0) { dx = 1; dy = 1; } // ~1.1 sigma out: near the rule's boundary + else if (n % 11 == 0) { dx = 3; dy = -2; } // disk edge: MINPK keeps it, the rule does not fire + else continue; + const int x = static_cast(std::lround(s.predicted[n].predicted_x)) + dx; + const int y = static_cast(std::lround(s.predicted[n].predicted_y)) + dy; + if (x < 0 || y < 0 || x >= static_cast(width) || y >= static_cast(height)) continue; + s.image[y * width + x] = (n % 2) ? INT32_MAX : INT32_MIN; + } return s; } @@ -126,7 +144,8 @@ void CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, float clip_nsigma = 4.0f, bool radial = false, int spacing = 60, float stencil_k = 0.0f, float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f, - OverlapMode overlap = OverlapMode::Off, float companion_dx = 0.0f) { + OverlapMode overlap = OverlapMode::Off, float companion_dx = 0.0f, + bool clip_spots = false) { const DiffractionExperiment experiment = MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial, DetJF(2), stencil_k, r1, r2, r3, overlap); @@ -135,7 +154,7 @@ void CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, const size_t npixel = experiment.GetPixelsNum(); REQUIRE(npixel == width * height); - const Scene scene = BuildScene(width, height, spacing, companion_dx); + const Scene scene = BuildScene(width, height, spacing, companion_dx, clip_spots); REQUIRE(scene.image.size() == npixel); REQUIRE(scene.predicted.size() > 60); @@ -161,6 +180,18 @@ void CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, // atomic summation of the learned profile, so compare up to a small tolerance. REQUIRE(out_gpu.size() == out_cpu.size()); REQUIRE(out_cpu.size() > 40); + if (clip_spots) { + // Guard against the coverage going vacuous: the punched pixels have to actually cost some + // reflections, or the two engines are being compared on a case neither of them meets. + const Scene clean_scene = BuildScene(width, height, spacing, companion_dx, false); + ImagePreprocessorBuffer clean_image(npixel); + for (size_t i = 0; i < npixel; ++i) + clean_image[i] = clean_scene.image[i]; + BraggIntegrationEngineCPU clean_cpu(experiment); + const auto out_clean = clean_cpu.Run(clean_image, clean_scene.predicted, + clean_scene.predicted.size(), 5); + CHECK(out_cpu.size() < out_clean.size()); + } for (size_t i = 0; i < out_cpu.size(); ++i) { INFO("mode " << static_cast(mode) << " reflection " << i << " hkl " << out_cpu[i].h); CHECK(out_gpu[i].h == out_cpu[i].h); @@ -239,6 +270,29 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") { 0.0f, 0.0f, 0.0f, OverlapMode::Exclude); } SECTION("ProfileGaussian mono trim") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 0.0f); } + // Unreadable pixels inside the signal disks themselves: the MINPK rescue keeps the reflection and + // fits it over what is left, and the peak-loss rule throws back the ones that lost the profile's + // maximum. Both decisions are per-reflection cuts on a reduction over the profile grid, computed + // independently in the two engines (serial max vs an atomicMax on the float bit pattern), so they + // have to reject exactly the same reflections - a mismatch shows up as a size mismatch here. + SECTION("ProfileGaussian clipped disks") { + CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 60, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, true); + } + SECTION("ProfileEmpirical clipped disks") { + CompareCpuVsGpu(IntegratorMode::ProfileEmpirical, std::nullopt, 4.0f, false, 60, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, true); + } + // The same, with an elongated profile: the peak is then a ridge, so the fraction-of-peak test has + // to protect a crest rather than one pixel, and the grid it reduces over is reflection-dependent. + SECTION("ProfileGaussian clipped disks stencil") { + CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.005f, 4.0f, false, 120, 3.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, true); + } + SECTION("BoxSum clipped disks") { + CompareCpuVsGpu(IntegratorMode::BoxSum, std::nullopt, 4.0f, false, 60, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, true); + } // The radial background curvature correction is computed independently in the two engines // (host loop vs radial_correct kernel), so it needs its own parity coverage. SECTION("BoxSum radial") { CompareCpuVsGpu(IntegratorMode::BoxSum, std::nullopt, 4.0f, true); }