diff --git a/tests/AdaptiveSpotFinderCPUTest.cpp b/tests/AdaptiveSpotFinderCPUTest.cpp index 7d8df5b7..f7819812 100644 --- a/tests/AdaptiveSpotFinderCPUTest.cpp +++ b/tests/AdaptiveSpotFinderCPUTest.cpp @@ -73,3 +73,58 @@ TEST_CASE("AdaptiveSpotFinderCPU_RawGeometry", "[AdaptiveSpotFinder]") { CHECK(std::lround(spots[0].RawCoord().x) == static_cast(spot_col)); CHECK(std::lround(spots[0].RawCoord().y) == static_cast(spot_row)); } + +// The property the whole engine exists for: the threshold comes from the image's OWN noise, so the +// same settings behave the same way on a frame whose background is ten times higher. A frame is built +// with background spread S around a mean, one pixel planted a few S above it (must stay unfound) and +// one planted far above (must be found); then the identical frame scaled by ten must give the identical +// answer. Any threshold that does not track the background - a constant, or one that drops the sigma +// term - finds the weak pixel in the scaled frame, or loses the strong one. +TEST_CASE("AdaptiveSpotFinderCPU_ThresholdTracksBackground", "[AdaptiveSpotFinder]") { + DiffractionExperiment x(DetJF4M()); + x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080); + x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.05, 5.0); + x.GeometryTransformation(false); + + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + const auto &pixel_to_bin = mapping.GetPixelToBin(); + + const size_t w = x.GetXPixelsNum(); + const size_t h = x.GetYPixelsNum(); + + // Two well-separated pixels that carry a ring, so both are seen by the finder. + std::vector planted; + for (size_t row = 300; row < h - 300 && planted.size() < 2; row += 137) + for (size_t col = 300; col < w - 300; col += 149) + if (pixel_to_bin[row * w + col] != UINT16_MAX) { + planted.push_back(row * w + col); + break; + } + REQUIRE(planted.size() == 2); + + // Background takes 5 evenly spaced levels one S apart, i.e. mean + 2S and sigma = sqrt(2) S. With + // ~100 expected noise pixels per frame the cut lands near mean + 4.1 sigma = mean + 5.8 S. + const auto run_at_scale = [&](int32_t scale) { + ImagePreprocessorBuffer buffer(x.GetPixelsNum()); + for (size_t i = 0; i < w * h; i++) + buffer[i] = scale * (10 + static_cast(i % 5)); + buffer[planted[0]] = scale * (10 + 5); // mean + 3 S: below the cut + buffer[planted[1]] = scale * (10 + 30); // mean + 28 S: well above it + std::vector res_mask(x.GetPixelsNum(), false); + AdaptiveSpotFinderCPU finder(mapping); + return finder.Run(buffer, AdaptiveSettings(), res_mask); + }; + + const auto plain = run_at_scale(1); + const auto scaled = run_at_scale(10); + + REQUIRE(plain.size() == 1); + CHECK(std::lround(plain[0].RawCoord().x) == static_cast(planted[1] % w)); + CHECK(std::lround(plain[0].RawCoord().y) == static_cast(planted[1] / w)); + + // Ten times the background, ten times the noise, ten times the signal - same answer. + REQUIRE(scaled.size() == plain.size()); + CHECK(std::lround(scaled[0].RawCoord().x) == std::lround(plain[0].RawCoord().x)); + CHECK(std::lround(scaled[0].RawCoord().y) == std::lround(plain[0].RawCoord().y)); +} diff --git a/tests/AdaptiveSpotFinderGPUTest.cpp b/tests/AdaptiveSpotFinderGPUTest.cpp index 2b86e994..ab01ebec 100644 --- a/tests/AdaptiveSpotFinderGPUTest.cpp +++ b/tests/AdaptiveSpotFinderGPUTest.cpp @@ -71,8 +71,8 @@ std::vector> SortedCoords(const std::vector } // namespace // Spot-finding functionality: the fused GPU engine must reproduce the reference CPU adaptive finder's -// spot list (the two share AdaptiveThreshold.h and the host connected-component extractor; the only -// difference is the GPU's float atomic ring reduction, which is exact for a realistic background). +// spot list. The two share AdaptiveThreshold.h and the host connected-component extractor, and both +// sum the rings in double, so the only difference left is the order the ring sums are accumulated in. TEST_CASE("AdaptiveSpotFinderGPU_SpotFindingParity", "[AdaptiveSpotFinderGPU]") { if (get_gpu_count() == 0) { WARN("No CUDA GPU present. Skipping AdaptiveSpotFinderGPU_SpotFindingParity"); @@ -147,6 +147,40 @@ TEST_CASE("AdaptiveSpotFinderGPU_AzimuthalIntegration", "[AdaptiveSpotFinderGPU] } } +// The ring sums are built by atomics, which arrive in an arbitrary order, so the same frame has to be +// re-run to show the engine agrees with itself: detection is a hard "value >= threshold" on integer +// counts, and a threshold that wobbles between runs flips pixels on the boundary and with them the size +// of a connected component. Two runs, same spot list. +TEST_CASE("AdaptiveSpotFinderGPU_RunToRunReproducible", "[AdaptiveSpotFinderGPU]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping AdaptiveSpotFinderGPU_RunToRunReproducible"); + return; + } + + DiffractionExperiment x = MakeExperiment(); + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + + ImagePreprocessorBufferGPU buffer(x.GetPixelsNum()); + FillTestImage(buffer, x); + REQUIRE(cudaMemcpy(buffer.getGPUBuffer(), buffer.getBuffer().data(), + x.GetPixelsNum() * sizeof(int32_t), cudaMemcpyHostToDevice) == cudaSuccess); + + std::vector res_mask(x.GetPixelsNum(), false); + const SpotFindingSettings settings = AdaptiveSettings(); + + auto stream = std::make_shared(); + AdaptiveSpotFinderGPU gpu(mapping, stream); + + const auto first = gpu.Run(buffer, settings, res_mask); + REQUIRE(first.size() > 0); + for (int repeat = 0; repeat < 4; repeat++) { + const auto again = gpu.Run(buffer, settings, res_mask); + REQUIRE(again.size() == first.size()); + REQUIRE(SortedCoords(again) == SortedCoords(first)); + } +} + TEST_CASE("AdaptiveSpotFinderGPU_Speed", "[AdaptiveSpotFinderGPU][.benchmark]") { if (get_gpu_count() == 0) { WARN("No CUDA GPU present. Skipping AdaptiveSpotFinderGPU_Speed");