jfjoch_viewer: Track top-pixels without storing/sorting all valid pixels

This commit is contained in:
2025-12-12 11:36:54 +01:00
parent d70c60c501
commit e03d0677ba
9 changed files with 241 additions and 51 deletions

View File

@@ -119,6 +119,8 @@ ADD_LIBRARY(JFJochCommon STATIC
ResolutionShells.h
DarkMaskSettings.cpp
DarkMaskSettings.h
TopPixels.cpp
TopPixels.h
)
TARGET_LINK_LIBRARIES(JFJochCommon JFJochLogger Compression JFCalibration gemmi Threads::Threads -lrt )

4
common/TopPixels.cpp Normal file
View File

@@ -0,0 +1,4 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "TopPixels.h"

86
common/TopPixels.h Normal file
View File

@@ -0,0 +1,86 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#ifndef JFJOCH_TOPPIXELS_H
#define JFJOCH_TOPPIXELS_H
#include <array>
#include <algorithm>
#include <cstdint>
#include <stdexcept>
#include <vector>
class TopPixels {
public:
struct Entry {
int32_t value = 0;
int32_t index = 0;
};
explicit TopPixels(int capacity = 20)
: capacity_(std::clamp(capacity, 1, kMaxCapacity)) {}
void Clear() noexcept { size_ = 0; }
void Add(int32_t value, int32_t index) noexcept {
const Entry e{value, index};
if (size_ < capacity_) {
InsertSorted_(e);
return;
}
// Quick reject: if not better than the current smallest in top-K
if (value <= top_[capacity_ - 1].value)
return;
ReplaceSmallestAndFixOrder_(e);
}
[[nodiscard]] int Capacity() const noexcept { return capacity_; }
[[nodiscard]] int Size() const noexcept { return size_; }
[[nodiscard]] bool Empty() const noexcept { return size_ == 0; }
// Sorted descending by value; valid for i in [0, Size()).
[[nodiscard]] const Entry& At(int i) const {
if (i < 0 || i >= size_)
throw std::out_of_range("TopPixels::At index out of range");
return top_[i];
}
[[nodiscard]] const Entry& operator[](int i) const { return At(i); }
// Convenience: copy results to a std::vector (still sorted descending).
[[nodiscard]] std::vector<Entry> ToVector() const {
return std::vector<Entry>(top_.begin(), top_.begin() + size_);
}
private:
static constexpr int kMaxCapacity = 64; // hard cap for stack-friendly storage
std::array<Entry, kMaxCapacity> top_{};
int size_ = 0;
int capacity_ = 20;
void InsertSorted_(Entry e) noexcept {
top_[size_] = e;
++size_;
// Bubble up to keep descending order by value
for (int i = size_ - 1; i > 0 && top_[i].value > top_[i - 1].value; --i)
std::swap(top_[i], top_[i - 1]);
}
void ReplaceSmallestAndFixOrder_(Entry e) noexcept {
// Replace the smallest (last, because sorted descending)
top_[capacity_ - 1] = e;
// Bubble up to restore descending order
for (int i = capacity_ - 1; i > 0 && top_[i].value > top_[i - 1].value; --i)
std::swap(top_[i], top_[i - 1]);
size_ = capacity_;
}
};
#endif //JFJOCH_TOPPIXELS_H

View File

@@ -6,6 +6,9 @@
#include "JFJochDecompress.h"
#include "../image_analysis/bragg_integration/BraggIntegrate2D.h"
#include <queue>
#include <algorithm>
JFJochReaderImage::JFJochReaderImage(const DataMessage &in_message,
const std::shared_ptr<const JFJochReaderDataset> &in_dataset)
: message(in_message),
@@ -72,11 +75,22 @@ void JFJochReaderImage::ProcessInputImage(const void *data, size_t npixel, int64
const T* img_ptr = reinterpret_cast<const T*>(data);
size_t good_pixel = 0;
// Reset per-image stats
saturated_pixel.clear();
error_pixel.clear();
valid_count = 0;
has_valid = false;
valid_pixel.reserve(image.size());
top_pixels_acc.Clear();
top_pixels.clear();
top_pixels.reserve(top_pixels_acc.Capacity());
for (int i = 0; i < npixel; i++) {
bool has_input_mask = false;
const auto &mask = dataset->pixel_mask.GetMask();
if (mask.size() == npixel)
has_input_mask = true;
for (size_t i = 0; i < npixel; i++) {
int32_t val;
if (img_ptr[i] <= INT32_MAX)
val = static_cast<int32_t>(img_ptr[i]);
@@ -84,44 +98,61 @@ void JFJochReaderImage::ProcessInputImage(const void *data, size_t npixel, int64
val = INT32_MAX;
uint32_t mask_val = 0;
if (!dataset->pixel_mask.GetMask().empty())
mask_val = dataset->pixel_mask.GetMask().at(i);
if (has_input_mask)
mask_val = mask[i];
if ((mask_val & (
(1<<PixelMask::ModuleGapPixelBit)
| (1<<PixelMask::ChipGapPixelBit)
| (1<<PixelMask::ModuleEdgePixelBit))) != 0) {
image[i] = GAP_PXL_VALUE;
} else if ((mask_val != 0) || (img_ptr[i] == special_value)){
} else if ((mask_val != 0) || (img_ptr[i] == special_value)) {
image[i] = ERROR_PXL_VALUE;
error_pixel.emplace(i);
error_pixel.emplace(static_cast<int64_t>(i));
} else if (val >= sat_value) {
image[i] = SATURATED_PXL_VALUE;
saturated_pixel.emplace(i);
saturated_pixel.emplace(static_cast<int64_t>(i));
} else {
good_pixel++;
image[i] = val;
valid_pixel.emplace_back(std::make_pair(val, i));
if (!has_valid) {
has_valid = true;
valid_min = val;
valid_max = val;
} else {
valid_min = std::min(valid_min, val);
valid_max = std::max(valid_max, val);
}
valid_count++;
top_pixels_acc.Add(val, static_cast<int32_t>(i));
}
}
// Sort based on the first element only
std::sort(valid_pixel.begin(), valid_pixel.end(),
[](const auto& a, const auto& b) {return a.first < b.first;});
CalcAutoContrast();
// For now: auto-foreground based purely on max valid element
auto_foreground = has_valid ? std::max(1, valid_max) : 10;
// Export top pixels (already sorted descending) into the existing vector interface
for (int i = 0; i < top_pixels_acc.Size(); i++) {
const auto &e = top_pixels_acc[i];
top_pixels.emplace_back(e.value, e.index);
}
}
void JFJochReaderImage::CalcAutoContrast() {
if (valid_pixel.empty())
auto_foreground = 10;
else {
auto it = valid_pixel.crbegin();
const auto index = static_cast<size_t>(valid_pixel.size() * auto_foreground_range);
std::advance(it, std::max<size_t>(1ULL, index) - 1);
// Now simplified: based on max element (valid_max)
auto_foreground = has_valid ? std::max(1, valid_max) : 10;
}
auto_foreground = std::max(1, it->first);
}
std::optional<std::pair<int32_t, int32_t>> JFJochReaderImage::ValidMinMax() const {
if (!has_valid)
return {};
return std::make_pair(valid_min, valid_max);
}
const std::vector<std::pair<int32_t, int32_t>> &JFJochReaderImage::GetTopPixels() const {
return top_pixels;
}
const DataMessage &JFJochReaderImage::ImageData() const {
@@ -144,10 +175,6 @@ const std::unordered_set<int64_t> &JFJochReaderImage::ErrorPixels() const {
return error_pixel;
}
const std::vector<std::pair<int32_t, int32_t>> &JFJochReaderImage::ValidPixels() const {
return valid_pixel;
}
const JFJochReaderDataset &JFJochReaderImage::Dataset() const {
if (!dataset)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
@@ -167,37 +194,56 @@ void JFJochReaderImage::AddImage(const JFJochReaderImage &other) {
error_pixel.clear();
saturated_pixel.clear();
valid_pixel.clear();
valid_pixel.reserve(image.size());
for (int i = 0; i < image.size(); i++) {
if (image[i] == GAP_PXL_VALUE || other.image[i] == GAP_PXL_VALUE)
valid_count = 0;
has_valid = false;
top_pixels_acc.Clear();
top_pixels.clear();
top_pixels.reserve(top_pixels_acc.Capacity());
for (size_t i = 0; i < image.size(); i++) {
if (image[i] == GAP_PXL_VALUE || other.image[i] == GAP_PXL_VALUE) {
image[i] = GAP_PXL_VALUE;
else if (image[i] == ERROR_PXL_VALUE || other.image[i] == ERROR_PXL_VALUE) {
} else if (image[i] == ERROR_PXL_VALUE || other.image[i] == ERROR_PXL_VALUE) {
image[i] = ERROR_PXL_VALUE;
error_pixel.emplace(i);
error_pixel.emplace(static_cast<int64_t>(i));
} else if (image[i] == SATURATED_PXL_VALUE || other.image[i] == SATURATED_PXL_VALUE) {
image[i] = SATURATED_PXL_VALUE;
saturated_pixel.emplace(i);
saturated_pixel.emplace(static_cast<int64_t>(i));
} else {
int64_t sum = static_cast<int64_t>(image[i]) + static_cast<int64_t>(other.image[i]);
if (sum <= INT32_MIN + 5) [[unlikely]] {
image[i] = ERROR_PXL_VALUE;
error_pixel.emplace(i);
error_pixel.emplace(static_cast<int64_t>(i));
} else if (sum > dataset->experiment.GetSaturationLimit()) [[unlikely]] {
image[i] = SATURATED_PXL_VALUE;
saturated_pixel.emplace(i);
saturated_pixel.emplace(static_cast<int64_t>(i));
} else {
image[i] = static_cast<int32_t>(sum);
valid_pixel.emplace_back(std::make_pair(sum, i));
const int32_t val = static_cast<int32_t>(sum);
image[i] = val;
if (!has_valid) {
has_valid = true;
valid_min = val;
valid_max = val;
} else {
valid_min = std::min(valid_min, val);
valid_max = std::max(valid_max, val);
}
valid_count++;
top_pixels_acc.Add(val, static_cast<int32_t>(i));
}
}
}
// Sort based on the first element only
std::sort(valid_pixel.begin(), valid_pixel.end(),
[](const auto& a, const auto& b) { return a.first < b.first; });
// For now: auto-foreground based purely on max valid element
auto_foreground = has_valid ? std::max(1, valid_max) : 10;
CalcAutoContrast();
for (int i = 0; i < top_pixels_acc.Size(); i++) {
const auto &e = top_pixels_acc[i];
top_pixels.emplace_back(e.value, e.index);
}
}
std::vector<float> JFJochReaderImage::GetAzInt1D() const {

View File

@@ -8,7 +8,9 @@
#include <unordered_set>
#include <map>
#include <mutex>
#include <optional>
#include "../common/TopPixels.h"
#include "JFJochReaderDataset.h"
#include "../common/CrystalLattice.h"
@@ -26,6 +28,16 @@ class JFJochReaderImage {
std::unordered_set<int64_t> error_pixel;
std::vector<std::pair<int32_t, int32_t>> valid_pixel;
// Fast stats without storing/sorting all valid pixels
int32_t valid_min = 0;
int32_t valid_max = 0;
size_t valid_count = 0;
bool has_valid = false;
// For overlay: track top pixels with a tiny O(K) structure; export to vector for UI
TopPixels top_pixels_acc{20};
std::vector<std::pair<int32_t, int32_t>> top_pixels;
constexpr static float auto_foreground_range = 0.001f;
int32_t auto_foreground;
void CalcAutoContrast();
@@ -43,9 +55,11 @@ public:
const std::vector<int32_t> &Image() const;
const std::unordered_set<int64_t> &SaturatedPixels() const;
const std::unordered_set<int64_t> &ErrorPixels() const;
const std::vector<std::pair<int32_t, int32_t>> &ValidPixels() const;
const JFJochReaderDataset &Dataset() const;
std::optional<std::pair<int32_t, int32_t>> ValidMinMax() const;
const std::vector<std::pair<int32_t, int32_t>> &GetTopPixels() const;
void AddImage(const JFJochReaderImage& other);
std::vector<float> GetAzInt1D() const;
std::vector<float> GetAzInt1D_BinToQ() const;

View File

@@ -62,6 +62,7 @@ ADD_EXECUTABLE(jfjoch_test
LatticeSearchTest.cpp
TimeTest.cpp
RotationIndexerTest.cpp
TopPixelsTest.cpp
)
target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJochReader JFJochWriter JFJochImageAnalysis JFJochCommon JFJochHLSSimulation JFJochPreview)

28
tests/TopPixelsTest.cpp Normal file
View File

@@ -0,0 +1,28 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include "../common/TopPixels.h"
TEST_CASE("TopPixels") {
TopPixels tp(5);
REQUIRE(tp.Empty());
REQUIRE(tp.Size() == 0);
for (int i = 50; i > 0; i--)
tp.Add(i*2, i);
tp.Add(500,245);
tp.Add(0,1221);
REQUIRE(tp.Size() == 5);
CHECK(tp[0].index == 245);
CHECK(tp[0].value == 500);
for (int i = 1; i < 4; i++) {
CHECK(tp[i].index == (50 - i + 1));
CHECK(tp[i].value == (50 - i + 1) * 2);
}
}

View File

@@ -290,11 +290,15 @@ void JFJochDiffractionImage::DrawBeamCenter() {
void JFJochDiffractionImage::DrawTopPixels() {
int i = 0;
for (auto iter = image->ValidPixels().crbegin();
iter != image->ValidPixels().rend() && i < show_highest_pixels;
iter++, i++)
DrawCross(iter->second % image->Dataset().experiment.GetXPixelsNum() + 0.5,
iter->second / image->Dataset().experiment.GetXPixelsNum() + 0.5, 15, 3);
for (const auto& p : image->GetTopPixels()) {
if (i >= show_highest_pixels)
break;
const int32_t idx = p.second;
DrawCross(idx % image->Dataset().experiment.GetXPixelsNum() + 0.5,
idx / image->Dataset().experiment.GetXPixelsNum() + 0.5, 15, 3);
i++;
}
}

View File

@@ -238,10 +238,15 @@ void JFJochViewerImageStatistics::loadImage(std::shared_ptr<const JFJochReaderIm
.arg(image->ImageData().spot_count_ice_rings.value_or(0))
.arg(image->ImageData().spot_count_indexed.value_or(0)));
text = QString("<b>%1</b> - <b>%2</b>")
.arg(image->ValidPixels().begin()->first)
.arg(image->ValidPixels().rbegin()->first);
valid_values->setText(text);
auto mm = image->ValidMinMax();
if (mm.has_value()) {
text = QString("<b>%1</b> - <b>%2</b>")
.arg(mm->first)
.arg(mm->second);
valid_values->setText(text);
} else {
valid_values->setText("N/A");
}
valid_values->setToolTip(QString("Error pixels: <b>%1</b><br/>Saturated pixels: <b>%2</b>")
.arg(image->ErrorPixels().size()).arg(image->SaturatedPixels().size()));