diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index bf18df563..c68e939b3 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -604,6 +604,12 @@ they are the strongest reflections the crystal gave, so they concentrate at low and recommends attenuating the beam or shortening the exposure (or a second, attenuated low-resolution sweep). +Stills are different: a still has no rocking event, and its overloaded observations are the fully +recorded ones of the strongest reflections, so they are kept in the merge as a lower bound - the box sum +with every saturated pixel at the saturation limit - and are not counted above. +**`OVERLOADED_REFLECTIONS_KEPT=`** counts them and **`OVERLOADED_REFLECTIONS_KEPT_PCT=`** gives their share +of the integrated reflections; the `OVERLOADS` warning fires from 0.05%. + **`OVERLAPPED_REFLECTIONS_PCT=`** is the share of the predicted reflections whose signal disk a nearer neighbour's prediction takes pixels from, as the integrator's overlap treatment (`--overlap`, default `exclude`) assigns them, and **`OVERLAPPED_REFLECTIONS_MAJOR_PCT=`** the share left with less than diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index ca860600b..d563684df 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -226,6 +226,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, probe_counts.profile_fallback += after.profile_fallback - before.profile_fallback; probe_counts.overlapped += after.overlapped - before.overlapped; probe_counts.overlapped_major += after.overlapped_major - before.overlapped_major; + probe_counts.overloads_kept += after.overloads_kept - before.overloads_kept; return ret; }; @@ -484,6 +485,7 @@ BraggIntegrationCounts MXAnalysisWithoutFPGA::BraggCounts() const { c.profile_fallback -= probe_counts.profile_fallback; c.overlapped -= probe_counts.overlapped; c.overlapped_major -= probe_counts.overlapped_major; + c.overloads_kept -= probe_counts.overloads_kept; return c; } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index bfcb4401f..0c7a5cd12 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -127,6 +127,8 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe if (overlap == OverlapMode::Exclude && mode == IntegratorMode::BoxSum) overlap = OverlapMode::Off; overlap_min_peak = settings.GetOverlapMinPeak(); + keep_overloads = !experiment.IsRotationIndexing(); + saturation_limit = static_cast(experiment.GetSaturationLimit()); claim = static_cast(R); inv_claim = 1.0f / claim; @@ -233,7 +235,7 @@ std::vector BraggIntegrationEngine::Finalize(const std::vector &panels, const std::vector &results, - int64_t image_number) const { + int64_t image_number) { // Sized to what is kept, not to what was predicted: this vector is the one every image holds to the // end of the pass, and a background ring lost to the neighbours drops a reflection here - on a // dense pattern more than half of them, which a reserve of npredicted would carry as dead capacity. @@ -258,6 +260,20 @@ std::vector BraggIntegrationEngine::Finalize(const std::vector(fr.saturated * (saturation_limit - fr.bkg)); + refl.sigma = std::sqrt(fr.var_bkg + std::max(refl.I, 1.0f)); + refl.overloaded = false; + ++counts.overloads_kept; + } refl.observed = true; if (polarization) refl.prescaling_corr /= geom.CalcAzIntPolarizationCorr(panels[i].plane_x, panels[i].plane_y, diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index bec850887..4c74b443d 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -141,6 +141,8 @@ struct BraggIntegrationCounts { // only (the crowding of the pattern); nothing reads them. uint64_t overlapped = 0; uint64_t overlapped_major = 0; + // Overloaded reflections of stills, kept in the merge as a lower bound (BraggIntegrationEngine::Finalize). + uint64_t overloads_kept = 0; }; // One reflection's extracted intensity, produced by the derived engine and turned into a @@ -159,6 +161,7 @@ struct BraggFitResult { bool ok = false; bool has_observed = false; bool overloaded = false; // a signal-disk pixel was saturated (Reflection::overloaded) + int saturated = 0; // signal-disk pixels at the saturation limit, left out of I }; class BraggIntegrationEngine { @@ -202,6 +205,11 @@ protected: float overlap_min_peak = 0.0f; // Reject: least clean profile fraction that is kept float claim = 0.0f, inv_claim = 0.0f; + // Overloaded reflections of stills are kept as a lower bound (Finalize), with every saturated pixel + // counted at the saturation limit, in the units of the preprocessed image. + bool keep_overloads = false; + double saturation_limit = 0.0; + // --- radial background curvature correction (BraggIntegrationSettings) --- // The disk and the annulus are concentric, so any background LINEAR in position cancels between // them; what survives is the curvature of the radial background. mean_annulus(B) - mean_disk(B) @@ -255,7 +263,7 @@ protected: std::vector Finalize(const std::vector &predicted, size_t npredicted, const std::vector &panels, const std::vector &results, - int64_t image_number) const; + int64_t image_number); public: explicit BraggIntegrationEngine(const DiffractionExperiment &experiment); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index e99c37f81..e3881dca6 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -281,6 +281,8 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (px == INT32_MAX || (px == INT32_MIN && on && !static_masked(static_cast(y) * W + x))) s.overloaded = 1; + if (px == INT32_MAX) + ++s.n_saturated; if (!valid(px)) continue; s.I_sum += px; s.I_sum_x += static_cast(x) * px; @@ -368,6 +370,7 @@ void BraggIntegrationEngineCPU::FinishPassA(PassA &a, size_t i, const Reflection out.r_bin = std::clamp(static_cast(std::lround(st.r0)), 0, a.n_rad - 1); out.k_bin = BraggStencilKernelIndex(st, n_kern); out.overloaded = s.overloaded != 0; + out.saturated = s.n_saturated; if (overlap != OverlapMode::Off && s.n_own < s.n_disk) { ++counts.overlapped; @@ -635,7 +638,7 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & if (overlap == OverlapMode::Reject && rh.n_own < overlap_min_peak * rh.n_disk) continue; results[i] = {static_cast(rh.I), static_cast(rh.sigma), static_cast(rh.bkg), static_cast(rh.obs_x), static_cast(rh.obs_y), - static_cast(rh.var_bkg), true, rh.has_obs, rh.overloaded}; + static_cast(rh.var_bkg), true, rh.has_obs, rh.overloaded, rh.saturated}; } return Finalize(predicted, npredicted, panels, results, image_number); } @@ -751,13 +754,14 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & fs = {}; 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. + // missing. It is passed on with its box sum and the flag, for a rotation merge to drop its + // rocking event (Reflection::overloaded) and for a still to be kept as a lower bound + // (Finalize); 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(rh.I), static_cast(rh.sigma), static_cast(rh.bkg), static_cast(rh.obs_x), static_cast(rh.obs_y), - static_cast(rh.var_bkg), true, false, true}; + static_cast(rh.var_bkg), true, false, true, rh.saturated}; continue; } const int sh = rh.shell < 0 ? 0 : rh.shell; @@ -866,7 +870,7 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & results[i] = {static_cast(I), static_cast(sigma), static_cast(rh.bkg), static_cast(rh.obs_x), static_cast(rh.obs_y), - static_cast(var_bkg), true, rh.has_obs, rh.overloaded}; + static_cast(var_bkg), true, rh.has_obs, rh.overloaded, rh.saturated}; } return Finalize(predicted, npredicted, panels, results, image_number); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h index 237fa4ad3..5e121a800 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h @@ -109,6 +109,7 @@ protected: 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 + int saturated = 0; // signal-disk pixels at the saturation limit }; // Everything pass A hands the profile passes. struct PassA { diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 2c1428552..7ba323174 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -166,7 +166,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p BraggGpuParams p, int n, int *cx_o, int *cy_o, float *I_o, float *sigma_o, float *bkg_o, float *bkgvar_o, float *varbkg_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o, - uint8_t *hasobs_o, uint8_t *overloaded_o, unsigned long long *invd2mm, + uint8_t *hasobs_o, uint8_t *overloaded_o, int *saturated_o, unsigned long long *invd2mm, float *isum_o, int *ninner_o, int *rbin_o, int *kbin_o, unsigned long long *rad_sum, int *rad_cnt, int n_rad, unsigned long long *counts) { @@ -175,7 +175,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p __shared__ unsigned long long s_Isum, s_Ix, s_Iy; __shared__ unsigned long long s_x, s_y; // positions behind s_Isum, to take the background out of the centroid - __shared__ int s_ninner, s_ninner_valid, s_nbkg, s_ndisk, s_nown, s_nsat; + __shared__ int s_ninner, s_ninner_valid, s_nbkg, s_ndisk, s_nown, s_nsat, s_nmax; __shared__ int s_nring, s_nclean; // readable ring pixels, and those no neighbour covers; see the CPU engine __shared__ unsigned long long s_bkgsum; __shared__ int s_accept, s_full; @@ -202,7 +202,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p s_r0 = sqrtf(rx * rx + ry * ry); s_Isum = 0; s_Ix = 0; s_Iy = 0; s_x = 0; s_y = 0; s_ninner = 0; s_ninner_valid = 0; s_nbkg = 0; s_bkgsum = 0; - s_ndisk = 0; s_nown = 0; s_nring = 0; s_nclean = 0; s_nsat = 0; + s_ndisk = 0; s_nown = 0; s_nring = 0; s_nclean = 0; s_nsat = 0; s_nmax = 0; } __syncthreads(); @@ -219,7 +219,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p long long l_Isum = 0, l_Ix = 0, l_Iy = 0; long long l_x = 0, l_y = 0; // positions behind l_Isum, for the background-free centroid - int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nring = 0, l_nclean = 0, l_nsat = 0; + int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nring = 0, l_nclean = 0, l_nsat = 0, l_nmax = 0; long long l_bkg = 0; for (int t = threadIdx.x; t < area; t += blockDim.x) { const int x = x0 + t % bw, y = y0 + t / bw; @@ -240,6 +240,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p const int32_t px = on ? img[y * p.W + x] : INT32_MIN; // Saturated, or unreadable on this frame beyond the run's mask: an overload (see the CPU engine). if (px == INT32_MAX || (px == INT32_MIN && on && !static_mask[y * p.W + x])) ++l_nsat; + if (px == INT32_MAX) ++l_nmax; if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px; l_x += x; l_y += y; ++l_niv; } } else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) { @@ -267,6 +268,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p WARP_ATOMIC_ADD(s_nring, l_nring); WARP_ATOMIC_ADD(s_nclean, l_nclean); WARP_ATOMIC_ADD(s_nsat, l_nsat); + WARP_ATOMIC_ADD(s_nmax, l_nmax); __syncthreads(); for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { s_radv[t] = 0; s_radn[t] = 0; } @@ -439,6 +441,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p kbin_o[i] = BraggStencilKernelIndex(st, p.n_kern); obsx_o[i] = (float) ox; obsy_o[i] = (float) oy; hasobs_o[i] = hasobs; overloaded_o[i] = s_nsat > 0 ? 1 : 0; + saturated_o[i] = s_nmax; ok_o[i] = 1; // The profile and its resolution shells are learned from COMPLETE reflections only; see the CPU // engine. @@ -945,6 +948,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { d_strong = CudaDevicePtr(new_capacity); d_has_obs = CudaDevicePtr(new_capacity); d_overloaded = CudaDevicePtr(new_capacity); + d_saturated = CudaDevicePtr(new_capacity); h_px_x = CudaHostPtr(new_capacity); h_px_y = CudaHostPtr(new_capacity); @@ -960,6 +964,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { h_ok = CudaHostPtr(new_capacity); h_has_obs = CudaHostPtr(new_capacity); h_overloaded = CudaHostPtr(new_capacity); + h_saturated = CudaHostPtr(new_capacity); capacity = new_capacity; } @@ -1048,7 +1053,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu boxsum<<>>(d_px_x, d_px_y, d_panel, d_d, img, d_mask, d_owner, d_static_mask->get(), p, n, d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y, - d_ok, d_strong, d_has_obs, d_overloaded, d_invd2, + d_ok, d_strong, d_has_obs, d_overloaded, d_saturated, d_invd2, d_isum, d_ninner, d_rbin, d_kbin, d_rad_sum, d_rad_cnt, rad_n, d_counts); cuda_err(cudaGetLastError()); @@ -1098,6 +1103,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu cuda_err(cudaMemcpyAsync(h_obs_y.get(), d_obs_y, sizeof(float) * npredicted, cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaMemcpyAsync(h_has_obs.get(), d_has_obs, sizeof(uint8_t) * npredicted, cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaMemcpyAsync(h_overloaded.get(), d_overloaded, sizeof(uint8_t) * npredicted, cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaMemcpyAsync(h_saturated.get(), d_saturated, sizeof(int) * npredicted, cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaStreamSynchronize(*stream)); // The clearing pass has run, so the buffers are clean again for the next image. if (clear_by_box) @@ -1111,6 +1117,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu results[i].var_bkg = h_var_bkg[i]; results[i].ok = true; results[i].overloaded = h_overloaded[i] != 0; + results[i].saturated = h_saturated[i]; if (h_has_obs[i]) { results[i].observed_x = h_obs_x[i]; results[i].observed_y = h_obs_y[i]; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index 8a4de8253..99ad82ab2 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -48,6 +48,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { // preprocessor's copy. See BraggIntegrationEngineCPU::static_mask. std::shared_ptr> d_static_mask; CudaDevicePtr d_ok, d_strong, d_has_obs, d_overloaded; + CudaDevicePtr d_saturated; // signal-disk pixels at the saturation limit // --- radial background curvature correction (see BraggIntegrationEngine) --- int n_rad = 0; // radial bins, 0 when the correction is off @@ -83,6 +84,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { CudaHostPtr h_mark; CudaHostPtr h_I, h_sigma, h_bkg, h_var_bkg, h_obs_x, h_obs_y; CudaHostPtr h_ok, h_has_obs, h_overloaded; + CudaHostPtr h_saturated; void EnsureCapacity(size_t n); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal index 46d12341f..d8f3316d0 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal +++ b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal @@ -125,6 +125,8 @@ kernel void box_sums(device const float2 *xy [[buffer(0)]], && ((static_mask[(uint) idx / 32] >> ((uint) idx % 32)) & 1u); if (v == PX_SATURATED || (v == PX_MASKED && on && !masked)) s.overloaded = 1; + if (v == PX_SATURATED) + ++s.n_saturated; if (!valid_px(v)) continue; s.I_sum += v; s.I_sum_x += (int64_t) x * v; diff --git a/image_analysis/bragg_integration/BraggPassA.h b/image_analysis/bragg_integration/BraggPassA.h index 8b1b10979..0193bab48 100644 --- a/image_analysis/bragg_integration/BraggPassA.h +++ b/image_analysis/bragg_integration/BraggPassA.h @@ -23,6 +23,8 @@ struct BraggBoxSums { int32_t n_bkg; // ring pixels behind bkg_sum int32_t n_ring, n_clean; // readable ring pixels, and those no predicted neighbour covers int32_t overloaded; // a signal-disk pixel was saturated, or unreadable on this frame only + int32_t n_saturated; // signal-disk pixels at the saturation limit + int32_t unused; }; // The ring pixels at or below a reflection's clip threshold (BraggIntegrationEngineCPU::ClipThreshold). diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index bc7bb3c64..a172fc49e 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1656,6 +1656,16 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, "collect a second, attenuated low-resolution sweep)", result.overloaded_reflections, pct, pct_low, OVERLOAD_LOW_RES_A)); } + if (result.overloads_kept > 0 && result.integrated_reflections > 0) { + const double pct = 100.0 * static_cast(result.overloads_kept) / result.integrated_reflections; + Add(s, KeyInt("OVERLOADED_REFLECTIONS_KEPT", static_cast(result.overloads_kept))); + Add(s, KeyReal("OVERLOADED_REFLECTIONS_KEPT_PCT", pct, "{:.3f}")); + if (pct >= OVERLOAD_WARNING_PCT) + Warn(doc, PathologyCode::OVERLOADS, fmt::format( + "Overloads: {} integrated reflections ({:.3f}%) held a saturated pixel. Stills keep them in the " + "merge as a lower bound, so the strongest reflections may read low. Attenuate the beam to " + "measure them", result.overloads_kept, pct)); + } if (result.predicted_reflections > 0 && experiment.GetBraggIntegrationSettings().GetOverlap() != OverlapMode::Off) { const double pct = 100.0 * static_cast(result.overlapped_reflections) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index b0d485b71..40ea445ca 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -787,6 +787,9 @@ struct ProcessResult { // predicted reflections whose signal disk a nearer neighbour takes pixels from, and of those the // ones left with less than --overlap-minpk of it. Zero with --overlap off. uint64_t predicted_reflections = 0, overlapped_reflections = 0, overlapped_reflections_major = 0; + // Overloaded reflections of stills, which are not flagged but kept as a lower bound + // (BraggIntegrationCounts::overloads_kept). + uint64_t overloads_kept = 0; // Which pass this result is (1-based) and how many ran, plus one line saying which pass became the // canonical output and on what evidence. pass_count == 1 outside the rotation two-pass. diff --git a/rugnux/RugnuxImagePass.cpp b/rugnux/RugnuxImagePass.cpp index 592bb2a77..de6d6442b 100644 --- a/rugnux/RugnuxImagePass.cpp +++ b/rugnux/RugnuxImagePass.cpp @@ -135,7 +135,8 @@ void Rugnux::ImagePass(PipelineLocals &p) { // finishes (BraggIntegrationCounts). Reading an engine's counts synchronises its stream, so it is // done once per worker at the end rather than per image. std::atomic bragg_predicted = 0, bragg_starved = 0, bragg_starved_nb = 0, - bragg_fallback = 0, bragg_overlapped = 0, bragg_overlapped_major = 0; + bragg_fallback = 0, bragg_overlapped = 0, bragg_overlapped_major = 0, + bragg_overloads_kept = 0; // Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and // the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of @@ -357,6 +358,7 @@ void Rugnux::ImagePass(PipelineLocals &p) { bragg_fallback += bragg.profile_fallback; bragg_overlapped += bragg.overlapped; bragg_overlapped_major += bragg.overlapped_major; + bragg_overloads_kept += bragg.overloads_kept; }; if (observer) { @@ -596,6 +598,11 @@ void Rugnux::ImagePass(PipelineLocals &p) { result.predicted_reflections = bragg_predicted; result.overlapped_reflections = bragg_overlapped; result.overlapped_reflections_major = bragg_overlapped_major; + result.overloads_kept = bragg_overloads_kept; + if (bragg_overloads_kept > 0) + logger.Info("Overloads: {} integrated reflections ({:.3f}% of the predicted) held a saturated pixel and are " + "kept as a lower bound, their saturated pixels at the saturation limit", + bragg_overloads_kept.load(), frac(bragg_overloads_kept)); if (experiment_.GetBraggIntegrationSettings().GetOverlap() != OverlapMode::Off) logger.Info("Overlapping spots: {:.2f}% of the predicted reflections share signal pixels with a " "nearer neighbour, {:.2f}% keep less than {:.0f}% of their profile as their own", diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index 067723b69..f5742d6d2 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -3,6 +3,7 @@ #include +#include #include #include #include @@ -112,12 +113,14 @@ TEST_CASE("BraggIntegrationEngineCPU_CrowdingCountedWithoutImage", "[Integration } } -// A saturated pixel (INT32_MAX, as the preprocessor marks it) at the peak of a spot: the reflection is -// not measured but kept, flagged overloaded, so that a rotation merge can drop its whole rocking event. -// The same spot without the saturated pixel is an ordinary measurement. +// A saturated pixel (INT32_MAX, as the preprocessor marks it) at the peak of a spot on rotation data: the +// reflection is not measured but kept, flagged overloaded, so that a rotation merge can drop its whole +// rocking event. The same spot without the saturated pixel is an ordinary measurement. TEST_CASE("BraggIntegrationEngineCPU_SaturatedPeakIsFlaggedNotDropped", "[Integration][portable]") { DiffractionExperiment experiment(DetJF(2)); - experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f); + experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f) + .Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {})); + experiment.ImportIndexingSettings(IndexingSettings().RotationIndexing(true)); experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings()); const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum(); @@ -148,3 +151,61 @@ TEST_CASE("BraggIntegrationEngineCPU_SaturatedPeakIsFlaggedNotDropped", "[Integr CHECK(saturated[0].overloaded); CHECK(saturated[0].clipped); } + +// A still has no rocking event to drop, so an overloaded observation is kept as a lower bound: the box sum +// with every saturated pixel at the saturation limit. A spot with its core clipped at the limit, the +// clipped pixels marked saturated as the preprocessor marks them: the still gives exactly the clipped +// sum, unflagged and clipped (it may be low), and below the unclipped spot; rotation data keep the flag. +TEST_CASE("BraggIntegrationEngineCPU_StillKeepsOverloadAsLowerBound", "[Integration][portable]") { + DiffractionExperiment experiment(DetJF(2)); + experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f); + experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings()); + const int32_t limit = 500, bkg = 12; + experiment.Detector().SaturationLimit(limit); + const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum(); + + const int cx = 600, cy = 300; + const float r1 = experiment.GetBraggIntegrationSettings().GetR1(); + ImagePreprocessorBuffer image(npixel); + for (size_t i = 0; i < npixel; ++i) + image[i] = bkg; + double clipped_sum = 0.0, true_sum = 0.0; + int n_saturated = 0; + for (int dy = -6; dy <= 6; ++dy) + for (int dx = -6; dx <= 6; ++dx) { + const int32_t v = bkg + static_cast(std::lround(2000.0 * std::exp(-(dx * dx + dy * dy) / (2 * 1.3 * 1.3)))); + image[(cy + dy) * width + cx + dx] = v >= limit ? INT32_MAX : v; + if (dx * dx + dy * dy < r1 * r1) { + clipped_sum += std::min(v, limit) - bkg; + true_sum += v - bkg; + n_saturated += v >= limit; + } + } + REQUIRE(n_saturated > 1); + + Reflection r{}; + r.h = 1; r.k = 2; r.l = 3; + r.predicted_x = cx; r.predicted_y = cy; + r.d = 2.0f; + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + const std::vector predicted{r}; + + BraggIntegrationEngineCPU still(experiment, PixelMask(experiment)); + const auto kept = still.Run(image, predicted, 1, 0); + REQUIRE(kept.size() == 1); + CHECK_FALSE(kept[0].overloaded); + CHECK(kept[0].clipped); + CHECK(kept[0].I == Catch::Approx(clipped_sum).epsilon(0.01)); + CHECK(kept[0].I < 0.8 * true_sum); + CHECK(kept[0].sigma >= std::sqrt(kept[0].I)); + CHECK(still.Counts().overloads_kept == 1); + + experiment.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {})) + .ImportIndexingSettings(IndexingSettings().RotationIndexing(true)); + BraggIntegrationEngineCPU rotation(experiment, PixelMask(experiment)); + const auto flagged = rotation.Run(image, predicted, 1, 0); + REQUIRE(flagged.size() == 1); + CHECK(flagged[0].overloaded); + CHECK(rotation.Counts().overloads_kept == 0); +} diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index c4030bc4f..f86d0af0a 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -204,9 +204,10 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, const auto out_clean = clean_cpu.Run(clean_image, clean_scene.predicted, clean_scene.predicted.size(), 5); // An unreadable pixel the (empty) mask does not explain reads as an overload, which keeps the - // reflection, flagged - so the cost shows in the measured ones. + // reflection - flagged, or on a still as a lower bound, clipped either way - so the cost shows + // in the complete ones. const auto measured = std::count_if(out_cpu.begin(), out_cpu.end(), - [](const Reflection &r) { return !r.overloaded; }); + [](const Reflection &r) { return !r.clipped; }); CHECK(static_cast(measured) < out_clean.size()); } for (size_t i = 0; i < out_cpu.size(); ++i) {