Azimuthal integration: clear sum2 along with sum and count
Clear() reallocated sum and count but left sum2 as the constructor had sized it. Every caller reuses one profile across the images of a dataset (rugnux keeps one per worker, the viewer one per view), so the squares kept accumulating while the means restarted at zero: the per-image standard deviation written to /entry/azint and shown in the plots was meaningless from the second image on, and grew without bound over a run. The size mismatch was the sharper edge. Clearing to a mapping with a different bin count left sum2 shorter than sum, and GetStd() and operator+= then read past its end - reachable in the viewer by opening a dataset with a wider q range than the one before it. Both are covered by tests: the same frame twice with a Clear() in between has to give the same standard deviation, and a profile cleared to a wider mapping has to report the right value in a bin that only the wider mapping has. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
@@ -44,6 +44,7 @@ void AzimuthalIntegrationProfile::Clear(const AzimuthalIntegrationMapping &mappi
|
||||
azim_bins = mapping.GetAzimuthalBinCount();
|
||||
|
||||
sum = std::vector<float>(mapping.GetBinNumber(), 0);
|
||||
sum2 = std::vector<float>(mapping.GetBinNumber(), 0);
|
||||
count = std::vector<uint64_t>(mapping.GetBinNumber(), 0);
|
||||
}
|
||||
|
||||
|
||||
@@ -256,6 +256,69 @@ TEST_CASE("AzimuthalIntegrationProfile_GetStd","[AzimuthalIntegration]") {
|
||||
CHECK(ret_count[3] == 0);
|
||||
}
|
||||
|
||||
TEST_CASE("AzimuthalIntegrationProfile_GetStd_AfterClear","[AzimuthalIntegration]") {
|
||||
DiffractionExperiment x(DetJF4M());
|
||||
x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
|
||||
x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);
|
||||
|
||||
PixelMask pixel_mask(x);
|
||||
AzimuthalIntegrationMapping mapping(x, pixel_mask);
|
||||
|
||||
AzimuthalIntegrationProfile profile(mapping);
|
||||
|
||||
std::vector<float> sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber());
|
||||
std::vector<uint32_t> count(mapping.GetBinNumber());
|
||||
|
||||
sum[0] = 2 + 3 + 4 + 5;
|
||||
sum2[0] = 4 + 9 + 16 + 25;
|
||||
count[0] = 4;
|
||||
|
||||
// The same frame twice, with a Clear() in between: a profile is reused for every image of a
|
||||
// dataset, so the second image has to give exactly the first one's standard deviation.
|
||||
profile.Add(sum, sum2, count);
|
||||
const auto first = profile.GetStd();
|
||||
|
||||
profile.Clear(mapping);
|
||||
profile.Add(sum, sum2, count);
|
||||
const auto second = profile.GetStd();
|
||||
|
||||
CHECK(first[0] == Catch::Approx(std::sqrt(5.0 / 3.0)));
|
||||
CHECK(second[0] == Catch::Approx(first[0]));
|
||||
}
|
||||
|
||||
TEST_CASE("AzimuthalIntegrationProfile_ClearToLargerMapping","[AzimuthalIntegration]") {
|
||||
DiffractionExperiment x(DetJF4M());
|
||||
x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
|
||||
x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4);
|
||||
|
||||
PixelMask pixel_mask(x);
|
||||
AzimuthalIntegrationMapping small(x, pixel_mask);
|
||||
|
||||
// A wider q range - the viewer re-uses one profile across datasets, so every vector Clear()
|
||||
// touches has to end up the size the new mapping asks for.
|
||||
x.QRangeForAzimInt_recipA(0.1, 9);
|
||||
AzimuthalIntegrationMapping large(x, pixel_mask);
|
||||
REQUIRE(large.GetBinNumber() > small.GetBinNumber());
|
||||
|
||||
AzimuthalIntegrationProfile profile(small);
|
||||
profile.Clear(large);
|
||||
|
||||
std::vector<float> sum(large.GetBinNumber()), sum2(large.GetBinNumber());
|
||||
std::vector<uint32_t> count(large.GetBinNumber());
|
||||
|
||||
// A bin that exists only in the wider mapping.
|
||||
const auto bin = small.GetBinNumber();
|
||||
sum[bin] = 2 + 3 + 4 + 5;
|
||||
sum2[bin] = 4 + 9 + 16 + 25;
|
||||
count[bin] = 4;
|
||||
|
||||
profile.Add(sum, sum2, count);
|
||||
|
||||
const auto stddev = profile.GetStd();
|
||||
REQUIRE(stddev.size() == large.GetBinNumber());
|
||||
CHECK(stddev[bin] == Catch::Approx(std::sqrt(5.0 / 3.0)));
|
||||
}
|
||||
|
||||
TEST_CASE("AzimuthalIntegrationProfile_operatorAdd","[AzimuthalIntegration]") {
|
||||
DiffractionExperiment x(DetJF4M());
|
||||
x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
|
||||
|
||||
Reference in New Issue
Block a user