diff --git a/common/Reflection.h b/common/Reflection.h index bcb614e40..a6b5ad9c5 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -58,6 +58,10 @@ struct Reflection { // estimate over what remained. Such a measurement may be low, and the merge does not let it // testify against a larger observation of the same reflection (WilsonOutliers.h). bool clipped = false; + // A pixel of the signal disk was saturated. The reflection is then not a measurement of its + // intensity at all - its brightest pixels are missing - and a rotation merge drops the whole + // rocking event it belongs to, as XDS does with an overloaded reflection. + bool overloaded = false; }; // One full reflection of a rotation sweep as the merge used it - its partials summed into one diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index a0229bed3..cb8a64565 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -9,6 +9,7 @@ * Rugnux: French-Wilson amplitudes (`F`/`SIGF`) use an anisotropic Wilson prior from the fitted anisotropy tensor, as ctruncate does; intensities are unchanged. * Rugnux claims a screw axis from a short axial row of a few weak reflections, so such crystals (e.g. P2_1 with a ~30 A unique axis) are no longer written without the screw. * Rugnux keeps the screw axes it found when a higher point group is adopted after the twin-law check, instead of writing the group without screws (e.g. P 4 2 2 for P4_2 2_1 2). +* Rugnux drops a rotation reflection whose spot held a saturated (overloaded) pixel, as XDS does, and reports the count as `OBSERVATIONS_REJECTED_OVERLOAD=`. ### 1.0.0-rc.173 diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index 17eda627b..66eb238c0 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -379,7 +379,8 @@ The combine groups each reflection's partials into rocking events (contiguous ru - **De-biased weighted sum.** Partials are combined by inverse-variance weighting, where each partial's variance is its background-noise component plus the *model* signal shared across the event (Kabsch profile-fit form). Using the shared model signal rather than the individual down-fluctuating intensity stops weak partials from being over-weighted, which would otherwise inflate the merged error model. The weights depend on the full, so the estimate is iterated. - **Captured fraction.** The partiality summed over the event, $f=\min(1,\sum_j p_j)$, measures how completely the rocking curve was sampled. A full whose curve was captured below a threshold (`--min-captured-fraction`, default 0.7 for rotation) is dropped — an event seen over only a small fraction of its curve is unreliable however many frames it spans. (The per-partial minimum-partiality cut of §10.2 still applies upstream, in the per-frame scaling.) - **Per-image rejection (opt-in).** A frame whose observations correlate poorly with the merged reference is not measuring the crystal being merged — it may be off-crystal, or on a *different* crystal where two lattices occupy separate regions of the sample. `--min-image-cc` drops such frames. It has no default: the per-frame correlation measures data quality as much as frame validity, and its typical level varies widely between datasets, so no single absolute bound is generally valid. -- **Capture-aware uncertainty.** A full captured incompletely ($f<1$) is extrapolated and biased high. The unobserved fraction is charged as an extra systematic uncertainty, $\sigma^2 \leftarrow \sigma^2 + \big(c\,(1-f)\,I\big)^2$, so the merge down-weights these extrapolated fulls and the error model treats their scatter as expected. It is enabled by default for the rotation path. +- **Capture-aware uncertainty.** A full captured incompletely ($f<1$) is extrapolated and biased high. The unobserved fraction is charged as an extra systematic uncertainty, $\sigma^2 \leftarrow \sigma^2 + \big(c\,(1-f)\,I\big)^2$, so the merge down-weights these extrapolated fulls and the error model treats their scatter as expected. The merge rebuilds every full's variance at the reflection's mean intensity (§10.4), and the capture term is rebuilt there too, as $\big(c\,(1-f)\,\langle I\rangle\big)^2$. It is enabled by default for the rotation path. +- **Overloaded events.** An event in which any partial had a saturated pixel in its signal disk — or a pixel unreadable on that frame alone, beyond the run's pixel mask, which is how a detector that writes its error value for a pixel it could not count reports an overload — is dropped whole, as XDS drops an overloaded reflection. The brightest part of such a rocking curve is exactly what is missing, so neither the sum of the remaining partials nor their extrapolation by the partiality model measures the reflection: on a strongly diffracting small-molecule crystal these were the strongest low-order reflections, and they read 2–3× low. The integration keeps an overloaded partial, unfitted and flagged, only so that the event can be recognised; nothing else reads it. The count is `OBSERVATIONS_REJECTED_OVERLOAD=` in the report. The fulls are then re-scaled in the XDS sense — a per-image scale refit directly on the complete reflections under the unity partiality model — and merged (§10.4). Because every merged observation is now a counting-statistics-limited full rather than a partiality-divided slice, the error model reaches a far higher asymptotic $I/\sigma$. diff --git a/image_analysis/MXAnalysisAfterFPGA.cpp b/image_analysis/MXAnalysisAfterFPGA.cpp index 5c8150b12..f92a009cd 100644 --- a/image_analysis/MXAnalysisAfterFPGA.cpp +++ b/image_analysis/MXAnalysisAfterFPGA.cpp @@ -36,7 +36,9 @@ MXAnalysisAfterFPGA::MXAnalysisAfterFPGA(const DiffractionExperiment &in_experim integration(in_integration), indexer(indexer), prediction(CreateBraggPrediction(experiment.IsRotationIndexing())), - bragg_engine(std::make_unique(in_experiment)) { + // No pixel mask reaches this path: the FPGA image marks its own masked pixels, so only a + // saturated pixel is read as an overload here. + bragg_engine(std::make_unique(in_experiment, PixelMask(in_experiment))) { if (experiment.IsSpotFindingEnabled()) find_spots = true; diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index ce3d2b173..6933bd473 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -53,7 +53,7 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp auto cpu_preprocessor = std::make_unique(in_experiment, in_mask); preprocessor_cpu = cpu_preprocessor.get(); preprocessor = std::move(cpu_preprocessor); - bragg_engine = std::make_unique(in_experiment); + bragg_engine = std::make_unique(in_experiment, in_mask); if (experiment.ROI().size() >= 1) roi = std::make_unique(experiment); #ifdef JFJOCH_USE_CUDA @@ -70,7 +70,7 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp // enable_fused_adaptive_gpu = true, so on the GPU path the copy is off in practice. preprocessor = std::make_unique(in_experiment, in_mask, stream, /*copy_image_to_host=*/!enable_fused_adaptive_gpu); - bragg_engine = std::make_unique(in_experiment, stream); + bragg_engine = std::make_unique(in_experiment, stream, in_mask); if (experiment.ROI().size() >= 1) roi = std::make_unique(experiment, stream); if (enable_fused_adaptive_gpu) { diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 281cf5629..013b74f15 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -733,8 +733,10 @@ std::vector> SumRockingEvents(const std::vector> SumRockingEvents(const std::vector BraggIntegrationEngine::Finalize(const std::vector(xpixel), static_cast(ypixel), 0), owner(static_cast(xpixel), static_cast(ypixel), BRAGG_OWNER_NONE) {} @@ -158,6 +160,14 @@ std::vector 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(npredicted); double inv_d2_min = std::numeric_limits::max(), inv_d2_max = 0.0; @@ -235,6 +245,9 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, else if (exclude) continue; } ++n_inner; + if (px == INT32_MAX + || (px == INT32_MIN && !static_masked(static_cast(y) * W + x))) + out.overloaded = true; if (!valid(px)) continue; I_sum += px; I_sum_x += static_cast(x) * px; @@ -269,7 +282,7 @@ std::vector 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 BraggIntegrationEngineCPU::RunImpl(const Sampler &img, 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}; + static_cast(rh.var_bkg), true, rh.has_obs, rh.overloaded}; } return Finalize(predicted, npredicted, results, image_number); } @@ -515,6 +528,16 @@ std::vector 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(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}; + continue; + } const int sh = rh.shell < 0 ? 0 : rh.shell; int Rf = R; @@ -668,7 +691,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, 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}; + static_cast(var_bkg), true, rh.has_obs, rh.overloaded}; } return Finalize(predicted, npredicted, results, image_number); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h index f125f391b..4060c2f06 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h @@ -4,6 +4,7 @@ #pragma once #include "BraggIntegrationEngine.h" +#include "../../common/PixelMask.h" class CompressedImage; @@ -52,13 +53,16 @@ class BraggIntegrationEngineCPU : public BraggIntegrationEngine { }; TiledFrame refl_mask; TiledFrame owner; + // The run's pixel mask, packed 32 pixels to a word (PixelMask::GetPackedMask). An unreadable pixel + // it does not explain was unreadable on this frame only - an overload (see Reflection::overloaded). + std::vector static_mask; template std::vector RunImpl(const Sampler &img, const std::vector &predicted, size_t npredicted, int64_t image_number); public: - explicit BraggIntegrationEngineCPU(const DiffractionExperiment &experiment); + BraggIntegrationEngineCPU(const DiffractionExperiment &experiment, const PixelMask &mask); using BraggIntegrationEngine::Run; // keep the preprocessed-buffer overload visible diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 3f4e3d4ca..2acdbd862 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include "BraggIntegrationEngineGPU.h" +#include "../indexing/CudaSharedTables.h" using namespace bragg_engine; @@ -156,11 +157,11 @@ __device__ __forceinline__ float profile_float(unsigned long long v) { // --- Pass A box-sum: rough I / background / centroid / strong flag, one block per reflection. --- __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, - const int32_t *img, const uint8_t *mask, const uint32_t *owner, + const int32_t *img, const uint8_t *mask, const uint32_t *owner, const uint8_t *static_mask, 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, unsigned long long *invd2mm, + uint8_t *hasobs_o, uint8_t *overloaded_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) { @@ -169,7 +170,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, __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; + __shared__ int s_ninner, s_ninner_valid, s_nbkg, s_ndisk, s_nown, s_nsat; __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; @@ -195,7 +196,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, 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_ndisk = 0; s_nown = 0; s_nring = 0; s_nclean = 0; s_nsat = 0; } __syncthreads(); @@ -212,7 +213,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, 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; + 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; 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; @@ -229,6 +230,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, } ++l_ni; const int32_t px = img[y * p.W + x]; + // Saturated, or unreadable on this frame beyond the run's mask: an overload (see the CPU engine). + if (px == INT32_MAX || (px == INT32_MIN && !static_mask[y * p.W + x])) ++l_nsat; 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) { @@ -254,6 +257,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, WARP_ATOMIC_ADD(s_nown, l_no); WARP_ATOMIC_ADD(s_nring, l_nring); WARP_ATOMIC_ADD(s_nclean, l_nclean); + WARP_ATOMIC_ADD(s_nsat, l_nsat); __syncthreads(); for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { s_radv[t] = 0; s_radn[t] = 0; } @@ -263,7 +267,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, // renormalises to the pixels it can read and Pass B cuts on how much of the profile survived // (XDS's MINPK). A box sum has no profile to renormalise with. See the CPU engine. s_full = (s_ninner_valid == s_ninner) ? 1 : 0; - s_accept = ((s_full || p.partial_ok) && s_nbkg > 5) ? 1 : 0; + // An overloaded reflection is kept whatever the mode, for its flag (see the CPU engine). + s_accept = ((s_full || p.partial_ok || s_nsat > 0) && s_nbkg > 5) ? 1 : 0; // A reflection the disk would have kept but the ring cannot support: dropped whole for want of // a background. And a ring every predicted neighbour would starve: the density the guard on a // widened radius reads. One atomic per such reflection, so nothing is paid where there is none. @@ -420,6 +425,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, rbin_o[i] = min(max((int) lroundf(s_r0), 0), n_rad > 0 ? n_rad - 1 : 0); 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; ok_o[i] = 1; // The profile and its resolution shells are learned from COMPLETE reflections only; see the CPU // engine. @@ -607,7 +613,7 @@ __global__ void radial_correct(const unsigned long long *rad_sum, const int *rad __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x, const float *px_y, const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, const float *I_seed, const float *sigma_seed, const float *bkg_a, - const float *bkgvar_a, const float *varbkg_seed, const uint8_t *ok_a, + const float *bkgvar_a, const float *varbkg_seed, const uint8_t *ok_a, const uint8_t *overloaded_a, const float *shell_P, const float *global_P, const float *sigma2_r, const float *sigma2_t, const int *shell_n, float *I_o, float *sigma_o, float *varbkg_o, uint8_t *ok_o, BraggGpuParams p, int n, @@ -623,6 +629,8 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x __shared__ int s_Rf, s_Gf; if (!ok_a[i]) { if (threadIdx.x == 0) ok_o[i] = 0; return; } + // An overloaded reflection keeps its box sum and is not fitted (see the CPU engine). + if (overloaded_a[i]) return; const int cx = cx_a[i], cy = cy_a[i]; const BraggStencil st = MakeBraggStencil(px_x[i], px_y[i], p.stencil); @@ -803,7 +811,8 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x } // namespace BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment &experiment, - std::shared_ptr stream) + std::shared_ptr stream, + const PixelMask &mask) : BraggIntegrationEngine(experiment), stream(std::move(stream)), d_mask(npixel), @@ -822,6 +831,8 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment // On the engine's stream (the member: the parameter was moved from), where the kernels counting // into it run - the NULL stream is not ordered before those. cuda_err(cudaMemsetAsync(d_counts, 0, sizeof(unsigned long long) * COUNT_SLOTS, *this->stream)); + d_static_mask = SharedDeviceTable(mask.GetBinaryMask().data(), npixel, mask.GetBinaryMask().data(), + mask.GetBinaryMaskChecksum(), *this->stream); // Fit profile grid: R for empirical / box, up to 3R (radially elongated) for the Gaussian. const int max_Rf = empirical ? R : 3 * R; @@ -899,6 +910,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { d_ok = CudaDevicePtr(new_capacity); d_strong = CudaDevicePtr(new_capacity); d_has_obs = CudaDevicePtr(new_capacity); + d_overloaded = CudaDevicePtr(new_capacity); h_px_x = CudaHostPtr(new_capacity); h_px_y = CudaHostPtr(new_capacity); @@ -912,6 +924,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { h_obs_y = CudaHostPtr(new_capacity); h_ok = CudaHostPtr(new_capacity); h_has_obs = CudaHostPtr(new_capacity); + h_overloaded = CudaHostPtr(new_capacity); capacity = new_capacity; } @@ -993,9 +1006,10 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu cuda_err(cudaGetLastError()); mark_mask<<>>(d_px_x, d_px_y, d_mark, d_mask, d_owner, p, n, MASK_FLUX); cuda_err(cudaGetLastError()); - boxsum<<>>(d_px_x, d_px_y, d_d, img, d_mask, d_owner, p, n, + boxsum<<>>(d_px_x, d_px_y, 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_invd2, + d_ok, d_strong, d_has_obs, d_overloaded, d_invd2, d_isum, d_ninner, d_rbin, d_kbin, d_rad_sum, d_rad_cnt, rad_n, d_counts); cuda_err(cudaGetLastError()); @@ -1020,7 +1034,8 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu d_shell_P, d_global_P, d_sigma2_r, d_sigma2_t, p); cuda_err(cudaGetLastError()); fit<<>>(img, d_owner, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, - d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_ok, d_shell_P, d_global_P, + d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_ok, d_overloaded, + d_shell_P, d_global_P, d_sigma2_r, d_sigma2_t, d_shell_n, d_I, d_sigma, d_var_bkg, d_ok, p, n, d_counts); cuda_err(cudaGetLastError()); @@ -1043,6 +1058,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu cuda_err(cudaMemcpyAsync(h_obs_x.get(), d_obs_x, sizeof(float) * npredicted, cudaMemcpyDeviceToHost, *stream)); 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(cudaStreamSynchronize(*stream)); // The clearing pass has run, so the buffers are clean again for the next image. if (clear_by_box) @@ -1055,6 +1071,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu results[i].bkg = h_bkg[i]; results[i].var_bkg = h_var_bkg[i]; results[i].ok = true; + results[i].overloaded = h_overloaded[i] != 0; 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 e6aa5b497..a86f968f2 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -9,6 +9,7 @@ #include "BraggIntegrationEngine.h" #include "../indexing/CUDAMemHelpers.h" +#include "../../common/PixelMask.h" // CUDA engine: reproduces BraggIntegrationEngineCPU up to floating-point precision. Each stage is a // kernel with one CUDA block per reflection cooperating over the small window via shared-memory @@ -42,7 +43,10 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { CudaDevicePtr d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y; CudaDevicePtr d_isum; // box-sum raw sum, for the radial correction CudaDevicePtr d_ninner, d_rbin, d_kbin; - CudaDevicePtr d_ok, d_strong, d_has_obs; + // The run's pixel mask, one byte per pixel (PixelMask::GetBinaryMask), shared with the + // preprocessor's copy. See BraggIntegrationEngineCPU::static_mask. + std::shared_ptr> d_static_mask; + CudaDevicePtr d_ok, d_strong, d_has_obs, d_overloaded; // --- radial background curvature correction (see BraggIntegrationEngine) --- int n_rad = 0; // radial bins, 0 when the correction is off @@ -76,12 +80,13 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { CudaHostPtr h_px_x, h_px_y, h_d; 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; + CudaHostPtr h_ok, h_has_obs, h_overloaded; void EnsureCapacity(size_t n); public: - BraggIntegrationEngineGPU(const DiffractionExperiment &experiment, std::shared_ptr stream); + BraggIntegrationEngineGPU(const DiffractionExperiment &experiment, std::shared_ptr stream, + const PixelMask &mask); std::vector Run(const ImagePreprocessorBuffer &image, const std::vector &predicted, size_t npredicted, int64_t image_number) override; diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index 5eaaa6b8c..299736af7 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -169,7 +169,7 @@ PostRefineObservations GatherPostRefineObservations(std::vector::max(), lmax = std::numeric_limits::min(); for (const auto &r : outcomes[o].reflections) - if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) { + if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f && !r.overloaded) { keep++; lmin = std::min(lmin, r.h); lmax = std::max(lmax, r.h); @@ -200,7 +200,7 @@ PostRefineObservations GatherPostRefineObservations(std::vector &h_count = hist[lo / chunk]; for (int o = lo; o < hi; o++) for (const auto &r : outcomes[o].reflections) - if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) + if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f && !r.overloaded) h_count[r.h - h_lo]++; }); @@ -227,7 +227,7 @@ PostRefineObservations GatherPostRefineObservations(std::vector fill = hist[lo / chunk]; for (int o = lo; o < hi; o++) { for (const auto &r : outcomes[o].reflections) { - if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f) continue; + if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f || r.overloaded) continue; const float ox = std::isfinite(r.observed_x) ? r.observed_x : NAN; const float oy = std::isfinite(r.observed_y) ? r.observed_y : NAN; pts[fill[r.h - h_lo]++] = Partial{static_cast(r.h), static_cast(r.k), diff --git a/image_analysis/scale_merge/HKLKey.cpp b/image_analysis/scale_merge/HKLKey.cpp index d2667e19b..3a08dc090 100644 --- a/image_analysis/scale_merge/HKLKey.cpp +++ b/image_analysis/scale_merge/HKLKey.cpp @@ -73,7 +73,7 @@ bool HKLKeyGenerator::IsSystematicallyAbsent(const Reflection &r) const { } bool AcceptReflection(const Reflection &r, std::optional d_min_limit, std::optional d_max_limit) { - if (!std::isfinite(r.I)) + if (!std::isfinite(r.I) || r.overloaded) // a saturated spot is not a measurement return false; if (!std::isfinite(r.d) || r.d <= 0.0f) return false; @@ -90,7 +90,7 @@ bool AcceptReflection(const Reflection &r, std::optional d_min_limit, st } bool AcceptReflection(const Reflection &r, double d_min_limit, double d_max_limit) { - if (!std::isfinite(r.I)) + if (!std::isfinite(r.I) || r.overloaded) // a saturated spot is not a measurement return false; if (!std::isfinite(r.d) || r.d <= 0.0f) return false; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 603141f08..6119ab318 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -185,6 +185,10 @@ struct MergeStatistics { // BETTER by every number it reports, so the count has to be visible or the failure is silent. // Zero when rejection is off. size_t n_observations_rejected = 0; + // Rocking events (rotation) left out of the merge because a pixel of their spot was saturated + // somewhere along the rocking curve - the measurement lacks its brightest part (XDS's OVERLOAD). + // Not part of n_observations_rejected, which counts what the outlier tests took. + size_t n_reflections_rejected_overloaded = 0; // The part of those the Wilson test removed (see WilsonOutliers.h), each listed, in the written range. std::vector wilson_rejected; double radiation_damage_delta_b = NAN; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 105f79550..ed6172a67 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -817,7 +817,7 @@ void RotationScaleMerge::BuildInRangeObservations(std::unique_ptr &ke fbuf[10][at] = r.predicted_y; frm[at] = o; onice[at] = r.on_ice_ring ? 1 : 0; - clip[at] = r.clipped ? 1 : 0; + clip[at] = (r.clipped ? OBS_CLIPPED : 0) | (r.overloaded ? OBS_OVERLOADED : 0); ++at; } }); @@ -854,7 +854,7 @@ void RotationScaleMerge::BuildInRangeObservations(std::unique_ptr &ke obs.image_number = r.image_number; obs.frame = o; obs.on_ice = r.on_ice_ring ? 1 : 0; - obs.clipped = r.clipped ? 1 : 0; + obs.clipped = (r.clipped ? OBS_CLIPPED : 0) | (r.overloaded ? OBS_OVERLOADED : 0); obs.corr = r.image_scale_corr; obs.group = -1; finite_ok[at] = (std::isfinite(obs.I) && std::isfinite(obs.prescaling_corr) && obs.prescaling_corr != 0.0f @@ -4149,11 +4149,12 @@ bool RotationScaleMerge::UsableFull(const Obs &o) const { void RotationScaleMerge::Combine() { fulls.clear(); g_full.assign(n_frames, 1.0); + fulls_dropped_overloaded = 0; // Combine one raw-hkl run into fulls, appended to `out`. Independent per run (the events of one hkl // touch no shared state), so runs parallelise cleanly. Returns the number of usable partials seen. // `dump` (serial path only) writes each emitted full for the diagnostic observation dump. - auto process_rawrun = [&](int r, std::vector &out, std::ofstream *dump) -> size_t { + auto process_rawrun = [&](int r, std::vector &out, std::ofstream *dump, int64_t &n_overloaded) -> size_t { const int lo = rawrun_start[r], hi = rawrun_start[r] + rawrun_count[r]; std::vector ev; // usable perm-indices of this raw hkl, in image-number order ev.reserve(hi - lo); @@ -4196,8 +4197,18 @@ void RotationScaleMerge::Combine() { float peak_px = partials[ev[i]].px, peak_py = partials[ev[i]].py; float peak_partiality = -1.0f; const bool on_ice = partials[ev[i]].on_ice; - bool clipped = false; - for (size_t m = i; m < kk; ++m) clipped = clipped || partials[ev[m]].clipped; + uint8_t flags = 0; + for (size_t m = i; m < kk; ++m) flags |= partials[ev[m]].clipped; + // A saturated pixel anywhere in the rocking curve means its brightest part was not measured: + // the event is not a measurement of the reflection, and summing what is left - or scaling + // the flanks up by the partiality model - reads it 2-3x low. Dropped whole, as XDS drops an + // overloaded reflection. + if (flags & OBS_OVERLOADED) { + ++n_overloaded; + i = kk; + continue; + } + const bool clipped = flags & OBS_CLIPPED; for (size_t m = i; m < kk; ++m) { const auto &r2 = partials[ev[m]]; const double sigma_corr = static_cast(r2.sigma) * r2.corr; @@ -4248,9 +4259,9 @@ void RotationScaleMerge::Combine() { continue; double sigma_full = 1.0 / std::sqrt(sum_w); + const double capture = capture_uncertainty_coeff * (1.0 - std::min(1.0, sum_partiality)); if (capture_uncertainty_coeff > 0.0) { - const double frac = std::min(1.0, sum_partiality); - const double extra = capture_uncertainty_coeff * (1.0 - frac) * std::max(0.0, F); + const double extra = capture * std::max(0.0, F); sigma_full = std::sqrt(sigma_full * sigma_full + extra * extra); } @@ -4260,6 +4271,7 @@ void RotationScaleMerge::Combine() { full.sigma = static_cast(sigma_full); full.var_bkg = static_cast(var_bkg_full); full.var_per_I = static_cast(var_bkg_full * var_bkg_full * sum_cwb); + full.capture = static_cast(capture); full.d = d; full.prescaling_corr = 1.0f; full.partiality = 1.0f; @@ -4294,7 +4306,7 @@ void RotationScaleMerge::Combine() { } fulls.reserve(n_run); for (int r = 0; r < n_run; ++r) - n_used += process_rawrun(r, fulls, dump.is_open() ? &dump : nullptr); + n_used += process_rawrun(r, fulls, dump.is_open() ? &dump : nullptr, fulls_dropped_overloaded); } else { // Parallel over contiguous rawrun chunks; concatenate the per-thread fulls in run order so the // result is deterministic. @@ -4302,6 +4314,7 @@ void RotationScaleMerge::Combine() { const int chunk = (n_run + nt - 1) / nt; std::vector> part(nt); std::vector used(nt, 0); + std::vector overloaded(nt, 0); std::vector> futures; futures.reserve(nt); for (int t = 0; t < nt; ++t) { @@ -4309,19 +4322,22 @@ void RotationScaleMerge::Combine() { if (r0 >= r1) break; futures.emplace_back(std::async(std::launch::async, [&, t, r0, r1] { part[t].reserve(r1 - r0); - for (int r = r0; r < r1; ++r) used[t] += process_rawrun(r, part[t], nullptr); + for (int r = r0; r < r1; ++r) used[t] += process_rawrun(r, part[t], nullptr, overloaded[t]); })); } for (auto &f : futures) f.get(); size_t total = 0; - for (int t = 0; t < nt; ++t) { total += part[t].size(); n_used += used[t]; } + for (int t = 0; t < nt; ++t) { + total += part[t].size(); n_used += used[t]; fulls_dropped_overloaded += overloaded[t]; + } fulls.reserve(total); for (int t = 0; t < nt; ++t) fulls.insert(fulls.end(), part[t].begin(), part[t].end()); } SortFullsByFrame(); - logger.Info("3D combine: {} fulls from {} partials", fulls.size(), n_used); + logger.Info("3D combine: {} fulls from {} partials, {} rocking events dropped for a saturated pixel", + fulls.size(), n_used, fulls_dropped_overloaded); } void RotationScaleMerge::SortFullsByFrame() { @@ -4525,10 +4541,14 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // below . Rebuild the variance at the reflection's EXPECTED intensity - var_bkg + var_per_I*, // the linear model the combine measured - so the weight no longer knows this full's own fluctuation. // Mirrors MergeOnTheFly::CorrectedSigma on the stills path. The error model is fitted on this same - // variance, so that its a scales what the merge applies it to. + // variance, so that its a scales what the merge applies it to. An event that caught only part of its + // rocking curve carries its capture uncertainty too (Obs::capture), at the same expected intensity: + // left out, a full extrapolated from three quarters of its curve merged at the weight of a whole one. auto counting_variance = [](const Obs &o, double I_for_b, double own_s2) { + const double I_exp = std::max(0.0, I_for_b); const double base = static_cast(o.corr) * o.corr * o.var_bkg - + static_cast(o.corr) * o.var_per_I * std::max(0.0, I_for_b); + + static_cast(o.corr) * o.var_per_I * I_exp + + static_cast(o.capture) * o.capture * I_exp * I_exp; return base > 0.0 ? base : own_s2; }; // Fit (a, b) to a pool of samples (ErrorModel.h). A lambda because the pool changes once the cutoff @@ -5690,6 +5710,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // Radiation-damage monitor (measured before any correction by MeasureRadiationDamageB): carry the // first->last relative-B change and the per-batch curve into the reported statistics / mmCIF. out.n_observations_rejected = reject_count; + out.n_reflections_rejected_overloaded = static_cast(fulls_dropped_overloaded); for (size_t j = 0; j < wilson.rejected.size(); ++j) { if (!wilson.rejected[j]) continue; const Obs &o = fulls[wilson_full[j]]; @@ -6108,7 +6129,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st if (gpu_active_ && observation_dump_path.empty()) { // The smoothed corr is already resident (scaling + smooth-G ran on the device, no round-trip). const int nf = gpu_->Combine(rawrun_group.data(), min_partiality, capture_uncertainty_coeff, - min_captured_fraction, max_frame_gap); + min_captured_fraction, max_frame_gap, fulls_dropped_overloaded); g_full.assign(n_frames, 1.0); if (scale_fulls && nf > 0) { @@ -6156,14 +6177,14 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st st.image_number.data(), st.frame.data(), st.on_ice.data(), st.clipped.data(), st.group.data()); gpu_->GetFullsPxPy(st.px.data(), st.py.data()); - gpu_->GetFullsVariance(st.var_bkg.data(), st.var_per_I.data()); + gpu_->GetFullsVariance(st.var_bkg.data(), st.var_per_I.data(), st.capture.data()); if (scaled_fulls_on_gpu) gpu_->GetFullsCorr(st.corr.data()); ParallelChunks(nf, ThreadsForWork(nf, nthreads), [&](int lo, int hi) { for (int i = lo; i < hi; ++i) { Obs &o = fulls[i]; o.h = st.h[i]; o.k = st.k[i]; o.l = st.l[i]; o.I = st.I[i]; o.sigma = st.sigma[i]; o.d = st.d[i]; - o.var_bkg = st.var_bkg[i]; o.var_per_I = st.var_per_I[i]; + o.var_bkg = st.var_bkg[i]; o.var_per_I = st.var_per_I[i]; o.capture = st.capture[i]; o.prescaling_corr = 1.0f; o.partiality = 1.0f; o.corr = scaled_fulls_on_gpu ? st.corr[i] : 1.0f; o.image_number = st.image_number[i]; o.frame = st.frame[i]; @@ -6174,7 +6195,8 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st o.zeta = 0.0f; o.delta_phi = 0.0f; o.bkg = 0.0f; } }); - logger.Info("3D combine{} (GPU): {} fulls", scaled_fulls_on_gpu ? " + scale-fulls" : "", nf); + logger.Info("3D combine{} (GPU): {} fulls, {} rocking events dropped for a saturated pixel", + scaled_fulls_on_gpu ? " + scale-fulls" : "", nf, fulls_dropped_overloaded); combined_on_gpu = true; } // The CPU combine is the one host stage left that reads corr off the partials, and it only runs for diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index a13e1be77..3baa4d15f 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -153,6 +153,9 @@ public: [[nodiscard]] double GetSearchMinZeta() const { return search_min_zeta; } private: + static constexpr uint8_t OBS_CLIPPED = 1; + static constexpr uint8_t OBS_OVERLOADED = 2; + // One integrated observation - a per-frame partial during scaling/combine, or a combined full during // scale-fulls/merge. Flat (not nested per image); a POD so the arrays translate straight to CUDA. struct Obs { @@ -161,11 +164,16 @@ private: // Fulls only, written by the combine: the full's variance as a function of intensity, // var(I) = var_bkg + var_per_I * I. The merge rebuilds it at the reflection's mean. float var_per_I = 0.0f; + // Fulls only: the capture uncertainty, a relative sigma - capture_uncertainty_coeff times the + // part of the rocking curve the event did not catch - so the merge can add (capture * )^2 to + // the variance it rebuilds at the reflection's mean, as the combine added it to the full's own. + float capture = 0.0f; float px = NAN, py = NAN; // predicted detector position (for the absorption surface; CPU path only) float image_number; // fractional frame position (for 3D-combine contiguity) int32_t frame; // index of the outcome whose per-frame scale G applies to this obs uint8_t on_ice; - uint8_t clipped; // Reflection::clipped; a full is clipped if any of its partials is + uint8_t clipped; // OBS_CLIPPED | OBS_OVERLOADED (Reflection::clipped / ::overloaded); a full + // is clipped if any of its partials is, and an overloaded event is no full float corr; // image_scale_corr (working; updated by scaling) int32_t group; // dense ASU-group id for the current space group; <0 = never mergeable }; @@ -345,6 +353,8 @@ private: void Resize(int n) { group.resize(n); I.resize(n); sigma.resize(n); corr.resize(n); d.resize(n); } }; MergeFields merge_fields; + // Rocking events the last combine dropped because a pixel of one of their partials was saturated. + int64_t fulls_dropped_overloaded = 0; // One host array per field for the fulls download: the device hands back an array per field and the // host gathers them into `fulls`. Members rather than locals in Run() because the whole @@ -352,12 +362,12 @@ private: // between them, so as locals every chain allocates, faults in and zeroes the lot again. struct FullsStaging { std::vector h, k, l, frame, group; - std::vector I, sigma, d, image_number, corr, px, py, var_bkg, var_per_I; + std::vector I, sigma, d, image_number, corr, px, py, var_bkg, var_per_I, capture; std::vector on_ice, clipped; void Resize(int n) { h.resize(n); k.resize(n); l.resize(n); frame.resize(n); group.resize(n); I.resize(n); sigma.resize(n); d.resize(n); image_number.resize(n); corr.resize(n); - px.resize(n); py.resize(n); var_bkg.resize(n); var_per_I.resize(n); + px.resize(n); py.resize(n); var_bkg.resize(n); var_per_I.resize(n); capture.resize(n); on_ice.resize(n); clipped.resize(n); } }; diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index 45b58a881..68536ff9b 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -271,6 +271,10 @@ namespace { return isfinite(I[i]) && isfinite(sigma[i]) && sigma[i] > 0.0f; } + // RotationScaleMerge::OBS_CLIPPED / OBS_OVERLOADED, the flag bits a partial's `clipped` byte carries. + constexpr uint8_t OBS_CLIPPED_FLAG = 1; + constexpr uint8_t OBS_OVERLOADED_FLAG = 2; + // All device pointers + scalars the combine kernels need, passed by value. struct CombineParams { int n_runs; @@ -287,9 +291,10 @@ namespace { *__restrict__ rr_h, *__restrict__ rr_k, *__restrict__ rr_l, *__restrict__ rr_group; int32_t *rr_nevents; // count pass output + int32_t *rr_noverloaded; // count pass output: events dropped for a saturated pixel const int32_t *rr_offset; // emit pass: per-run base offset into the fulls arrays int32_t *f_h, *f_k, *f_l, *f_frame, *f_group; - float *f_I, *f_sigma, *f_d, *f_img, *f_px, *f_py, *f_var_bkg, *f_var_per_I; + float *f_I, *f_sigma, *f_d, *f_img, *f_px, *f_py, *f_var_bkg, *f_var_per_I, *f_capture; uint8_t *f_on_ice, *f_clipped; }; @@ -304,7 +309,7 @@ namespace { const int hi = lo + p.rr_count[r]; const int group = p.rr_group[r]; - int n_emit = 0; + int n_emit = 0, n_overloaded = 0; int cursor = lo; while (cursor < hi) { while (cursor < hi && !CombineUsable(p.perm[cursor], p.I, p.sigma, p.corr)) ++cursor; @@ -351,8 +356,11 @@ namespace { float peak_px = p.px[first], peak_py = p.py[first]; float peak_partiality = -1.0f; const bool on_ice = p.on_ice[first]; - bool clipped = false; - for (int m = ev_start; m <= ev_end; ++m) clipped = clipped || p.clipped[p.perm[m]]; + uint8_t flags = 0; + for (int m = ev_start; m <= ev_end; ++m) flags |= p.clipped[p.perm[m]]; + // An overloaded event is not a measurement of the reflection; see the host combine. + if (flags & OBS_OVERLOADED_FLAG) { ++n_overloaded; continue; } + const bool clipped = flags & OBS_CLIPPED_FLAG; for (int m = ev_start; m <= ev_end; ++m) { const int i = p.perm[m]; if (!CombineUsable(i, p.I, p.sigma, p.corr)) continue; @@ -406,9 +414,9 @@ namespace { continue; double sigma_full = 1.0 / sqrt(sum_w); + const double capture = p.capture_uncertainty_coeff * (1.0 - Dmin(1.0, sum_partiality)); if (p.capture_uncertainty_coeff > 0.0) { - const double frac = Dmin(1.0, sum_partiality); - const double extra = p.capture_uncertainty_coeff * (1.0 - frac) * Dmax(0.0, F); + const double extra = capture * Dmax(0.0, F); sigma_full = sqrt(sigma_full * sigma_full + extra * extra); } @@ -419,6 +427,7 @@ namespace { p.f_sigma[o] = float(sigma_full); p.f_var_bkg[o] = float(var_bkg_full); p.f_var_per_I[o] = float(var_bkg_full * var_bkg_full * sum_cwb); + p.f_capture[o] = float(capture); p.f_d[o] = d; p.f_img[o] = peak_frame; p.f_px[o] = peak_px; p.f_py[o] = peak_py; @@ -430,8 +439,10 @@ namespace { ++n_emit; } - if (!Emit) + if (!Emit) { p.rr_nevents[r] = n_emit; + p.rr_noverloaded[r] = n_overloaded; + } } template @@ -451,7 +462,7 @@ namespace { int n_groups; double min_partiality, error_model_a, error_model_b, reject_nsigma; int for_search, error_model_active, reject_outliers; - const float *I, *sigma, *corr, *partiality, *d, *reject_median, *var_bkg, *var_per_I; + const float *I, *sigma, *corr, *partiality, *d, *reject_median, *var_bkg, *var_per_I, *capture; const float *reject_var_add; // per group: the shell's measured share of dI^2/4 const int32_t *group, *frame; const uint8_t *on_ice, *frame_cell_ok, *half; @@ -554,7 +565,9 @@ namespace { const double bi = p.error_model_b * I_for_b; const double c = p.corr[i]; double a_var = double(sigma_raw) * sigma_raw; - const double base = c * c * p.var_bkg[i] + c * p.var_per_I[i] * Dmax(0.0, I_for_b); + const double I_exp = Dmax(0.0, I_for_b); + const double cap = p.capture[i]; + const double base = c * c * p.var_bkg[i] + c * p.var_per_I[i] * I_exp + cap * cap * I_exp * I_exp; if (base > 0.0) a_var = base; const double v = p.error_model_a * a_var + bi * bi; return v > 0.0 ? float(sqrt(v)) : sigma_raw; @@ -900,11 +913,11 @@ struct RotationScaleMergeGPU::Impl { CudaDevicePtr bkg, var_bkg, image_number, d_obs, px_obs, py_obs; int n_runs = 0, n_perm = 0; CudaDevicePtr perm, rr_start, rr_count, rr_h, rr_k, rr_l, rr_group; - CudaDevicePtr rr_nevents, rr_offset; + CudaDevicePtr rr_nevents, rr_noverloaded, rr_offset; // combine: resident fulls SoA (rebuilt each Combine) int n_fulls = 0; CudaDevicePtr f_h, f_k, f_l, f_frame, f_group; - CudaDevicePtr f_I, f_sigma, f_d, f_img, f_px, f_py, f_var_bkg, f_var_per_I; + CudaDevicePtr f_I, f_sigma, f_d, f_img, f_px, f_py, f_var_bkg, f_var_per_I, f_capture; CudaDevicePtr f_on_ice, f_clipped; // scale-fulls (Unity model, kept resident): all-ones partiality/prescaling_corr/zeta so the shared scaling kernels // yield coeff=mean, plus the working corr, the per-obs scale scratch, and the fulls frame/group CSRs @@ -1215,7 +1228,7 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.half = d.m_half.get(); - p.var_bkg = d.f_var_bkg.get(); p.var_per_I = d.f_var_per_I.get(); + p.var_bkg = d.f_var_bkg.get(); p.var_per_I = d.f_var_per_I.get(); p.capture = d.f_capture.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); p.em_mean = d.m_em_mean.get(); p.cc_factor = d.cc_factor.get(); p.a_swI = d.a_swI.get(); p.a_sw = d.a_sw.get(); p.a_swIh0 = d.a_swIh0.get(); p.a_swIh1 = d.a_swIh1.get(); @@ -1270,6 +1283,7 @@ void RotationScaleMergeGPU::MergeRmeas(const double *merged_I, double *absdev, d p.r_sumv = d.r_sumv.get(); p.r_sumv2 = d.r_sumv2.get(); p.error_model_a = d.merge_em_a; p.error_model_b = d.merge_em_b; p.error_model_active = d.merge_em_active; p.em_mean = d.m_em_mean.get(); p.var_bkg = d.f_var_bkg.get(); p.var_per_I = d.f_var_per_I.get(); + p.capture = d.f_capture.get(); p.rejected_obs = d.m_rejected.get(); const int grp_blocks = std::min(65535, (ng + BLK - 1) / BLK); @@ -1373,12 +1387,13 @@ void RotationScaleMergeGPU::SetRawRuns(int n_runs, int n_perm, const int32_t *pe d.rr_group = d.Alloc(std::max(1, n_runs)); d.rr_nevents = d.Alloc(std::max(1, n_runs)); + d.rr_noverloaded = d.Alloc(std::max(1, n_runs)); d.rr_offset = d.Alloc(std::max(1, n_runs)); } int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_partiality, double capture_uncertainty_coeff, double min_captured_fraction, - float max_frame_gap) { + float max_frame_gap, int64_t &n_overloaded) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; CopyAndWait(d.rr_group.get(), rawrun_group, size_t(d.n_runs) * sizeof(int32_t), @@ -1397,6 +1412,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti p.perm = d.perm.get(); p.rr_start = d.rr_start.get(); p.rr_count = d.rr_count.get(); p.rr_h = d.rr_h.get(); p.rr_k = d.rr_k.get(); p.rr_l = d.rr_l.get(); p.rr_group = d.rr_group.get(); p.rr_nevents = d.rr_nevents.get(); + p.rr_noverloaded = d.rr_noverloaded.get(); const int blocks = std::min(65535, (d.n_runs + BLK - 1) / BLK); @@ -1412,6 +1428,11 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti int64_t acc = 0; for (int r = 0; r < d.n_runs; ++r) { offset[r] = static_cast(acc); acc += nevents[r]; } d.n_fulls = static_cast(acc); + std::vector noverloaded(d.n_runs); + CopyAndWait(noverloaded.data(), d.rr_noverloaded.get(), size_t(d.n_runs) * sizeof(int32_t), + cudaMemcpyDeviceToHost, impl_->s(), "download noverloaded"); + n_overloaded = 0; + for (int r = 0; r < d.n_runs; ++r) n_overloaded += noverloaded[r]; // Allocate the fulls SoA and emit. const int nf = std::max(1, d.n_fulls); @@ -1420,7 +1441,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti d.f_I = d.Alloc(nf); d.f_sigma = d.Alloc(nf); d.f_d = d.Alloc(nf); d.f_img = d.Alloc(nf); d.f_px = d.Alloc(nf); d.f_py = d.Alloc(nf); - d.f_var_bkg = d.Alloc(nf); d.f_var_per_I = d.Alloc(nf); + d.f_var_bkg = d.Alloc(nf); d.f_var_per_I = d.Alloc(nf); d.f_capture = d.Alloc(nf); d.f_on_ice = d.Alloc(nf); d.f_clipped = d.Alloc(nf); d.f_corr = d.Alloc(nf); d.f_partiality = d.Alloc(nf); d.f_rlp = d.Alloc(nf); d.f_zeta = d.Alloc(nf); @@ -1434,7 +1455,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti p.f_frame = d.f_frame.get(); p.f_group = d.f_group.get(); p.f_I = d.f_I.get(); p.f_sigma = d.f_sigma.get(); p.f_d = d.f_d.get(); p.f_img = d.f_img.get(); p.f_px = d.f_px.get(); p.f_py = d.f_py.get(); - p.f_var_bkg = d.f_var_bkg.get(); p.f_var_per_I = d.f_var_per_I.get(); + p.f_var_bkg = d.f_var_bkg.get(); p.f_var_per_I = d.f_var_per_I.get(); p.f_capture = d.f_capture.get(); p.f_on_ice = d.f_on_ice.get(); p.f_clipped = d.f_clipped.get(); if (d.n_fulls > 0) { CombineKernel<<s()>>>(p); @@ -1555,7 +1576,7 @@ void RotationScaleMergeGPU::GetFullsPxPy(float *px, float *py) const { CopyAndWait(py, d.f_py.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_py"); } -void RotationScaleMergeGPU::GetFullsVariance(float *var_bkg, float *var_per_I) const { +void RotationScaleMergeGPU::GetFullsVariance(float *var_bkg, float *var_per_I, float *capture) const { DeviceGuard guard(impl_->device, impl_->available); const auto &d = *impl_; if (d.n_fulls == 0) return; @@ -1563,6 +1584,7 @@ void RotationScaleMergeGPU::GetFullsVariance(float *var_bkg, float *var_per_I) c CopyAndWait(var_bkg, d.f_var_bkg.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_var_bkg"); CopyAndWait(var_per_I, d.f_var_per_I.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_var_per_I"); + CopyAndWait(capture, d.f_capture.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_capture"); } void RotationScaleMergeGPU::SetFullsCorr(const float *corr) { diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index 1fd87d5df..4a75a7874 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -143,9 +143,10 @@ public: // adds the capture-uncertainty term. rawrun_group (length n_runs) is the current space group's ASU // id per raw hkl (it becomes the full's group). Deterministic: fulls are emitted in raw-run-major, // event order (a count pass -> host prefix sum -> emit-at-offset), matching the CPU path. Returns the - // number of fulls (call GetFulls with buffers of that length). + // number of fulls (call GetFulls with buffers of that length); n_overloaded receives how many rocking + // events were dropped for a saturated pixel (see the host combine). int Combine(const int32_t *rawrun_group, double min_partiality, double capture_uncertainty_coeff, - double min_captured_fraction, float max_frame_gap); + double min_captured_fraction, float max_frame_gap, int64_t &n_overloaded); // Download the combined fulls SoA (length = Combine()'s return). The working corr is downloaded // separately by GetFullsCorr (it is only meaningful after the fulls scaling; otherwise the caller sets it). @@ -157,8 +158,8 @@ public: // surface. Length = n_fulls. void GetFullsPxPy(float *px, float *py) const; - // Download the fulls' variance model, var(I) = var_bkg + var_per_I * I. Length = n_fulls. - void GetFullsVariance(float *var_bkg, float *var_per_I) const; + // Download the fulls' variance model, var(I) = var_bkg + var_per_I * I + (capture * I)^2. Length = n_fulls. + void GetFullsVariance(float *var_bkg, float *var_per_I, float *capture) const; // Re-upload the fulls' working corr (length n_fulls) after the host correction surfaces (decay / // absorption) mutate it, so the resident merge reads the corrected scale. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 14e9848f4..7bfa100cf 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -902,6 +902,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, // Outlier rejection drops observations from the merge AND from R_meas and the CC1/2 // half-sets, so a run that rejects too much scores BETTER on every other number here. Add(s, KeyInt("OBSERVATIONS_REJECTED", ms.n_observations_rejected)); + // Rocking events whose spot held a saturated pixel: not a measurement, left out (XDS's OVERLOAD). + Add(s, KeyInt("OBSERVATIONS_REJECTED_OVERLOAD", ms.n_reflections_rejected_overloaded)); // The part of those removed as improbable under Wilson statistics (an observation with no // symmetry mates to be judged against, or out-voting them); each is listed in the developer // report. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 596ee0b0e..c7adb42a3 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -6766,7 +6766,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b for (const auto &io : indexer->GetIntegrationOutcome()) for (const auto &r : io.reflections) if ((r.h == 0 || r.k == 0 || r.l == 0) && std::isfinite(r.I) && std::isfinite(r.sigma) - && r.sigma > 0.0f) { + && r.sigma > 0.0f && !r.overloaded) { auto &sum = summed[{r.h, r.k, r.l}]; sum.first += r.I; sum.second += static_cast(r.sigma) * r.sigma; diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index f00690f42..8c68bb5e6 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -48,7 +48,7 @@ TEST_CASE("BraggIntegrationEngineCPU_NeighbourMaskFollowsPartiality", "[Integrat r.partiality = (gx == 0 && gy == 0) ? 0.5f : neighbour_partiality; predicted.push_back(r); } - BraggIntegrationEngineCPU engine(experiment); + BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment)); const auto out = engine.Run(image, predicted, predicted.size(), 0); counts = engine.Counts(); std::vector spot; @@ -76,3 +76,40 @@ TEST_CASE("BraggIntegrationEngineCPU_NeighbourMaskFollowsPartiality", "[Integrat CHECK(on_frame.bkg_starved_by_neighbour == on_frame.bkg_starved); CHECK(on_frame.bkg_starved_by_neighbour == tails.bkg_starved_by_neighbour); } + +// 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. +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.ImportBraggIntegrationSettings(BraggIntegrationSettings()); + const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum(); + + const int cx = 600, cy = 300; + ImagePreprocessorBuffer image(npixel); + for (size_t i = 0; i < npixel; ++i) + image[i] = 12; + for (int dy = -6; dy <= 6; ++dy) + for (int dx = -6; dx <= 6; ++dx) + image[(cy + dy) * width + cx + dx] += static_cast(std::lround(800.0 * std::exp(-(dx * dx + dy * dy) / (2 * 1.3 * 1.3)))); + + 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 engine(experiment, PixelMask(experiment)); + const auto clean = engine.Run(image, predicted, 1, 0); + REQUIRE(clean.size() == 1); + CHECK_FALSE(clean[0].overloaded); + + image[cy * width + cx] = INT32_MAX; + const auto saturated = engine.Run(image, predicted, 1, 0); + REQUIRE(saturated.size() == 1); + CHECK(saturated[0].overloaded); + CHECK(saturated[0].clipped); +} diff --git a/tests/BraggIntegrationEngineCompressedImageTest.cpp b/tests/BraggIntegrationEngineCompressedImageTest.cpp index 881c4b4b8..584924187 100644 --- a/tests/BraggIntegrationEngineCompressedImageTest.cpp +++ b/tests/BraggIntegrationEngineCompressedImageTest.cpp @@ -101,7 +101,7 @@ TEST_CASE("BraggIntegrationEngineCPU_CompressedImageMatchesBuffer", "[Integratio const Scene scene = BuildScene(W, H); REQUIRE(scene.image.size() == npixel); - BraggIntegrationEngineCPU engine(experiment); + BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment)); // Route A: the preprocessed int32 buffer. ImagePreprocessorBuffer buffer(npixel); diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index 55d160c6e..c4030bc4f 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -2,6 +2,8 @@ // SPDX-License-Identifier: GPL-3.0-only #include + +#include #include "../common/CUDAWrapper.h" #ifdef JFJOCH_USE_CUDA @@ -171,7 +173,7 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, ImagePreprocessorBuffer cpu_image(npixel); for (size_t i = 0; i < npixel; ++i) cpu_image[i] = scene.image[i]; - BraggIntegrationEngineCPU cpu(experiment); + BraggIntegrationEngineCPU cpu(experiment, PixelMask(experiment)); const auto out_cpu = cpu.Run(cpu_image, scene.predicted, scene.predicted.size(), 5); // GPU under test, identical input uploaded to the device @@ -181,7 +183,7 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, gpu_image[i] = scene.image[i]; REQUIRE(cudaMemcpyAsync(gpu_image.getGPUBuffer(), gpu_image.getBuffer().data(), npixel * sizeof(int32_t), cudaMemcpyHostToDevice, *stream) == cudaSuccess); - BraggIntegrationEngineGPU gpu(experiment, stream); + BraggIntegrationEngineGPU gpu(experiment, stream, PixelMask(experiment)); const auto out_gpu = gpu.Run(gpu_image, scene.predicted, scene.predicted.size(), 5); // The ok/observed decisions are deterministic geometry, so both engines return the same set in @@ -198,15 +200,20 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, ImagePreprocessorBuffer clean_image(npixel); for (size_t i = 0; i < npixel; ++i) clean_image[i] = clean_scene.image[i]; - BraggIntegrationEngineCPU clean_cpu(experiment); + BraggIntegrationEngineCPU clean_cpu(experiment, PixelMask(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()); + // 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. + const auto measured = std::count_if(out_cpu.begin(), out_cpu.end(), + [](const Reflection &r) { return !r.overloaded; }); + CHECK(static_cast(measured) < 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); CHECK(out_gpu[i].image_number == out_cpu[i].image_number); + CHECK(out_gpu[i].overloaded == out_cpu[i].overloaded); CHECK(out_gpu[i].bkg == Catch::Approx(out_cpu[i].bkg).epsilon(0.02).margin(0.5)); CHECK(out_gpu[i].I == Catch::Approx(out_cpu[i].I).epsilon(0.03).margin(2.0)); CHECK(out_gpu[i].sigma == Catch::Approx(out_cpu[i].sigma).epsilon(0.03).margin(0.5)); @@ -386,11 +393,11 @@ TEST_CASE("BraggIntegrationEngineGPU_ReusedEngineMatchesFresh") { }; auto stream_fresh = std::make_shared(); - BraggIntegrationEngineGPU fresh(experiment, stream_fresh); + BraggIntegrationEngineGPU fresh(experiment, stream_fresh, PixelMask(experiment)); const auto out_fresh = integrate(fresh, second, stream_fresh); auto stream_reused = std::make_shared(); - BraggIntegrationEngineGPU reused(experiment, stream_reused); + BraggIntegrationEngineGPU reused(experiment, stream_reused, PixelMask(experiment)); integrate(reused, first, stream_reused); const auto out_reused = integrate(reused, second, stream_reused); @@ -425,7 +432,7 @@ TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") { REQUIRE(npixel == width * height); auto stream = std::make_shared(); - BraggIntegrationEngineGPU gpu(experiment, stream); + BraggIntegrationEngineGPU gpu(experiment, stream, PixelMask(experiment)); for (int spacing : {28, 60}) { const Scene scene = BuildScene(width, height, spacing); const size_t nrefl = scene.predicted.size(); @@ -446,7 +453,7 @@ TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") { const auto t1 = std::chrono::steady_clock::now(); const double ms = std::chrono::duration(t1 - t0).count() / iters; - BraggIntegrationEngineCPU cpu(experiment); + BraggIntegrationEngineCPU cpu(experiment, PixelMask(experiment)); ImagePreprocessorBuffer cpu_image(npixel); for (size_t i = 0; i < npixel; ++i) cpu_image[i] = scene.image[i]; const auto c0 = std::chrono::steady_clock::now(); diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 408bce534..6dcb3c6dc 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -117,6 +117,11 @@ TEST_CASE("AcceptReflection_ResolutionLimits") { CHECK(AcceptReflection(r, 0.0, 0.0)); CHECK_FALSE(AcceptReflection(r, 0.0, 15.0)); CHECK(AcceptReflection(r, 2.0, 50.0)); + + // A reflection with a saturated pixel is no measurement, whatever its resolution. + r.overloaded = true; + CHECK_FALSE(AcceptReflection(r, std::nullopt, std::nullopt)); + CHECK_FALSE(AcceptReflection(r, 0.0, 0.0)); } // --- Completeness denominator --------------------------------------------------------------------