diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index 76b0752ae..71494dc0b 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -177,6 +177,17 @@ float BraggIntegrationSettings::GetBackgroundPixelGain() const { return bkg_pixel_gain; } +BraggIntegrationSettings &BraggIntegrationSettings::PixelGain(float input) { + check_finite("Pixel gain", input); + check_min("Pixel gain", input, 1.0); + pixel_gain = input; + return *this; +} + +float BraggIntegrationSettings::GetPixelGain() const { + return pixel_gain; +} + BraggIntegrationSettings &BraggIntegrationSettings::BackgroundRadialCorrection(std::optional input) { bkg_radial_correction = input; return *this; diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 9c8023add..18425d85b 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -103,6 +103,13 @@ class BraggIntegrationSettings { // there the photon clip cuts the background's upper tail at about 2 sigma and adds a positive offset // to every reflection. Set by rugnux from the pre-scan for SMV and marCCD input (BackgroundGain.h). float bkg_pixel_gain = 1.0f; + // Detector units per photon of the integrated intensities. Above 1 the variances the integrator + // forms - the box sum's, the profile fit's and the background's - are taken at this many ADU per + // photon, so the strong-spot selection and the profile/summation guard compare in photons, and + // each reflection leaves the integrator converted to photons (I, sigma, bkg, var_bkg), which is + // what the merge's own variance model (var_bkg + I) assumes. Set by rugnux for rotation sweeps of a + // CCD from the same measurement as bkg_pixel_gain; 1 for a detector that counts photons. + float pixel_gain = 1.0f; // Radial background curvature correction. The signal disk and the background annulus are // concentric, so for ANY background linear in position their means are equal - a plane fit buys // nothing and the leading error is the CURVATURE of the radial background, which the flat annulus @@ -181,6 +188,7 @@ public: BraggIntegrationSettings& BackgroundTrimFraction(float input); BraggIntegrationSettings& BackgroundClipNSigma(float input); BraggIntegrationSettings& BackgroundPixelGain(float input); + BraggIntegrationSettings& PixelGain(float input); BraggIntegrationSettings& BackgroundRadialCorrection(std::optional input); BraggIntegrationSettings& MaxHKL(std::optional input); BraggIntegrationSettings& Overlap(OverlapMode input); @@ -202,6 +210,7 @@ public: [[nodiscard]] float GetBackgroundTrimFraction() const; [[nodiscard]] float GetBackgroundClipNSigma() const; [[nodiscard]] float GetBackgroundPixelGain() const; + [[nodiscard]] float GetPixelGain() const; // Unset = auto (gate per image on the smooth-ice score); see bkg_radial_correction. [[nodiscard]] std::optional GetBackgroundRadialCorrection() const; [[nodiscard]] std::optional GetMaxHKL() const; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index 25769d19f..417febd8c 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -118,6 +118,7 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe // longer overrides it for broadband data. bkg_clip_nsigma = settings.GetBackgroundClipNSigma(); bkg_pixel_gain = settings.GetBackgroundPixelGain(); + pixel_gain = settings.GetPixelGain(); bkg_trim = bkg_clip_nsigma > 0.0f ? 0.0f : settings.GetBackgroundTrimFraction(); // Overlap treatment. Ownership is decided out to the fit grid's half size, which is where the @@ -265,10 +266,11 @@ std::vector BraggIntegrationEngine::Finalize(const std::vector px] float bkg_clip_nsigma; // high-outlier background sigma-clip multiplier (0 = no clip) float bkg_pixel_gain = 1.0f; // background pixel variance per unit mean the clip scales by + // ADU per photon (BraggIntegrationSettings::pixel_gain): every variance the engine forms is taken + // at it, and Finalize converts the reflections to photons with it. 1 changes nothing, bit for bit. + float pixel_gain = 1.0f; bool use_ellipse; // radially elongate the per-reflection Gaussian double c_radial; // radial variance coefficient of tan^2(2theta): parallax + capture diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index a3b432758..11a1459e0 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -428,9 +428,13 @@ void BraggIntegrationEngineCPU::FinishPassA(PassA &a, size_t i, const Reflection // error enters n_inner times over: var(I) = I_sum + n_inner^2 * bkg/n_bkg_used. Leaving // the second term out understates sigma by sqrt(1 + n_inner/n_bkg) - 1.109x at the // default r1=4/r2=6/r3=10 stencil, on every reflection of every dataset. + // bkg_var is in counted units, as the profile fit reads it; every variance formed from it is + // taken at the gain, so a CCD's sigma - and the strong flag and the guard read from it - is in + // photons rather than in ADU counted as photons. + const double g = pixel_gain; out.bkg_var = out.bkg / n_bkg_used; - out.var_bkg = static_cast(s.n_inner_valid) * out.bkg - + static_cast(s.n_inner_valid) * s.n_inner_valid * out.bkg_var; + out.var_bkg = g * (static_cast(s.n_inner_valid) * out.bkg + + static_cast(s.n_inner_valid) * s.n_inner_valid * out.bkg_var); out.I_sum = s.I_sum; out.n_inner = static_cast(s.n_inner_valid); out.n_disk = s.n_disk; @@ -438,7 +442,7 @@ void BraggIntegrationEngineCPU::FinishPassA(PassA &a, size_t i, const Reflection const double var_bkg_term = static_cast(s.n_inner_valid) * s.n_inner_valid * out.bkg_var; out.sigma = 1.0; if (s.I_sum > 0) { - out.sigma = std::max(out.sigma, std::sqrt(static_cast(s.I_sum) + var_bkg_term)); + out.sigma = std::max(out.sigma, std::sqrt(g * (static_cast(s.I_sum) + var_bkg_term))); // The centroid of the SIGNAL, not of the disk. Weighting by the raw counts weights by // signal plus background, and the background is flat over a disk centred on the // PREDICTION - so its own centroid is the prediction exactly, and it pulls the answer @@ -850,15 +854,17 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & // 1/den is the profile-fit variance with the background taken as exact. The fit is // I = sum(P*(px-bkg)/v) / sum(P^2/v), so dI/dbkg = -wsum/den and the background estimate's // own error adds (wsum/den)^2 * var(bkg) - the same term the box sum was missing. - double sigma = std::sqrt(1.0 / den + (wsum / den) * (wsum / den) * rh.bkg_var); + // Both variances at the gain (FinishPassA), so the guard below compares like with like. + const double g = pixel_gain; + double sigma = std::sqrt(g * (1.0 / den + (wsum / den) * (wsum / den) * rh.bkg_var)); // var_bkg is the NON-SIGNAL part of that variance, and 1/den is the fit's variance taken at the // fitted intensity, so the signal part to remove is I itself - not max(0, I). Clamping it leaves // a down-fluctuated reflection's own (deflated) variance standing as its background variance, // which is 2-3x too small; the merge then weights exactly the down-fluctuated observations up. // The whole point of a separate var_bkg is that it does not move with the observation's own // fluctuation, and 1/den - I is what does not (1/den grows with I one for one). - double var_bkg = std::max(0.0, 1.0 / den - I - + (wsum / den) * (wsum / den) * rh.bkg_var); + double var_bkg = g * std::max(0.0, 1.0 / den - I + + (wsum / den) * (wsum / den) * rh.bkg_var); // The seed is a sum over the disk the box sum actually read, so when pixels are missing from // both - excluded to a neighbour, or unreadable - the fit's full-profile intensity has to be // scaled down to that same disk before the two are comparable. Nothing dropped gives 1. diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index a99faa4d9..c802cf73b 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -20,6 +20,7 @@ struct BraggGpuParams { int R, G, GG; float bkg_clip_nsigma; // high-side background sigma-clip multiplier (0 = no clip) float bkg_pixel_gain; // background pixel variance per unit mean (1 = counted photons) + float pixel_gain; // ADU per photon every variance is taken at (see the CPU engine) int empirical; // ProfileEmpirical vs ProfileGaussian int use_ellipse; float bw_sigma; @@ -421,7 +422,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p double sigma = 1.0; uint8_t hasobs = 0; double ox = 0.0, oy = 0.0; if (Isum > 0) { - sigma = fmax(sigma, sqrt((double) Isum + var_bkg_term)); + sigma = fmax(sigma, sqrt((double) p.pixel_gain * ((double) Isum + var_bkg_term))); // The centroid of the SIGNAL. Weighting by raw counts weights by signal plus background, // and the background is flat over a disk centred on the PREDICTION, so its centroid is the // prediction and it pulls the answer there - the displacement comes out shrunk by @@ -441,7 +442,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p cy_o[i] = (int) lroundf(cy); I_o[i] = (float) I; sigma_o[i] = (float) sigma; bkg_o[i] = (float) bkg; bkgvar_o[i] = (float) bkg_var; - varbkg_o[i] = (float) ((double) s_ninner_valid * bkg + var_bkg_term); + varbkg_o[i] = (float) ((double) p.pixel_gain * ((double) s_ninner_valid * bkg + var_bkg_term)); isum_o[i] = (float) Isum; ninner_o[i] = s_ninner_valid; rbin_o[i] = min(max((int) lroundf(s_r0), 0), n_rad > 0 ? n_rad - 1 : 0); @@ -828,9 +829,10 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x // 1/s_den takes the background as exact; dI/dbkg = -s_wsum/s_den adds the background // estimate's own error (see the CPU engine). const float wr = s_wsum / s_den; - float I = s_I, sigma = sqrtf(1.0f / s_den + wr * wr * bkgvar_a[i]); + // Both at the gain, as the seed's (see the CPU engine). + float I = s_I, sigma = sqrtf(p.pixel_gain * (1.0f / s_den + wr * wr * bkgvar_a[i])); // The signal part to remove is I itself, not max(0, I) - see the CPU engine. - float var_bkg = fmaxf(0.0f, 1.0f / s_den - I + wr * wr * bkgvar_a[i]); + float var_bkg = p.pixel_gain * fmaxf(0.0f, 1.0f / s_den - I + wr * wr * bkgvar_a[i]); // Guard against profile-fit runaways (see the CPU engine): fall back to the summation seed // when the profile result diverges from it. Pixels missing from both - excluded to a // neighbour, or unreadable - scale the full-profile intensity back to the disk the seed @@ -1012,6 +1014,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu .R = R, .G = G, .GG = GG, .bkg_clip_nsigma = mode != IntegratorMode::BoxSum ? bkg_clip_nsigma : 0.0f, .bkg_pixel_gain = bkg_pixel_gain, + .pixel_gain = pixel_gain, .empirical = empirical ? 1 : 0, .use_ellipse = use_ellipse ? 1 : 0, .bw_sigma = static_cast(bw_sigma), .c_radial = static_cast(c_radial), diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 6247f61bb..97f259287 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -606,6 +606,10 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru if (g) { BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); bis.BackgroundPixelGain(std::max(1.0f, *g)); + // On a rotation sweep the same gain sets the integrated intensities' variances (and their + // unit, photons). Stills keep their sigmas in ADU for now. + if (experiment_.IsRotationIndexing()) + bis.PixelGain(std::max(1.0f, *g)); experiment_.ImportBraggIntegrationSettings(bis); } if (!g) @@ -614,8 +618,10 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru else logger.Info("Background gain: a background pixel scatters {:.2f}x its mean in ADU ({} frames){}", *g, (ordinals.size() + gain_frame_stride - 1) / gain_frame_stride, - *g > 1.0f ? fmt::format(" => background clip at mean + {:.1f}*sqrt({:.2f}*mean)", - experiment_.GetBraggIntegrationSettings().GetBackgroundClipNSigma(), *g) + *g > 1.0f ? fmt::format(" => background clip at mean + {:.1f}*sqrt({:.2f}*mean){}", + experiment_.GetBraggIntegrationSettings().GetBackgroundClipNSigma(), *g, + experiment_.IsRotationIndexing() + ? "; intensities and sigmas in photons at this gain" : "") : std::string(" - at or below 1, the photon clip is kept")); } diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index e4789cc3e..153ddb7bb 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -351,3 +351,46 @@ TEST_CASE("BraggIntegrationEngineCPU_StillKeepsOnlySaturatedOverloads", "[Integr CHECK(out[1].overloaded); CHECK(still.Counts().overloads_kept == 1); } + +// A detector that reports gain ADU per photon: the same spot read at gain 4 leaves the integrator in +// photons - a quarter of the ADU intensity, at the Poisson sigma of that many photons, i.e. half the +// sigma the ADU counted as photons gave. +TEST_CASE("BraggIntegrationEngineCPU_PixelGainGivesPhotons", "[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); + 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] = 40; + 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}; + + for (const auto mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian}) { + BraggIntegrationSettings settings; + settings.Integrator(mode); + experiment.ImportBraggIntegrationSettings(settings); + BraggIntegrationEngineCPU adu(experiment, PixelMask(experiment)); + const auto out_adu = adu.Run(image, predicted, 1, 0); + settings.PixelGain(4.0f); + experiment.ImportBraggIntegrationSettings(settings); + BraggIntegrationEngineCPU photons(experiment, PixelMask(experiment)); + const auto out_ph = photons.Run(image, predicted, 1, 0); + REQUIRE(out_adu.size() == 1); + REQUIRE(out_ph.size() == 1); + CHECK(out_ph[0].I == Catch::Approx(out_adu[0].I / 4.0f)); + CHECK(out_ph[0].sigma == Catch::Approx(out_adu[0].sigma / 2.0f)); + CHECK(out_ph[0].bkg == Catch::Approx(out_adu[0].bkg / 4.0f)); + CHECK(out_ph[0].var_bkg == Catch::Approx(out_adu[0].var_bkg / 4.0f)); + } +} diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index 233449cc5..de125599f 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -150,13 +150,15 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f, OverlapMode overlap = OverlapMode::Off, float companion_dx = 0.0f, bool clip_spots = false, float odd_partiality = 1.0f, - const SpotFootprint *footprint = nullptr) { + const SpotFootprint *footprint = nullptr, float pixel_gain = 1.0f) { DiffractionExperiment experiment = MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial, DetJF(2), stencil_k, r1, r2, r3, overlap); - if (footprint) { + { BraggIntegrationSettings settings = experiment.GetBraggIntegrationSettings(); - settings.Footprint(*footprint); + if (footprint) + settings.Footprint(*footprint); + settings.PixelGain(pixel_gain); experiment.ImportBraggIntegrationSettings(settings); } const size_t width = experiment.GetXPixelsNum(); @@ -324,6 +326,16 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 8, 0.0f, 0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 0.01f); } + // A CCD's gain: every variance at it, so the strong flag and the guard decide in photons, and the + // reflections leave in photons - computed independently in the two engines. + SECTION("ProfileGaussian pixel gain") { + CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 60, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 1.0f, nullptr, 4.0f); + } + SECTION("BoxSum pixel gain") { + CompareCpuVsGpu(IntegratorMode::BoxSum, std::nullopt, 4.0f, false, 60, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 1.0f, nullptr, 4.0f); + } 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