mx/Jungfraujoch
0
0
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 37 Wiki Activity

File writer and spot finding improvements

Browse Source
This commit is contained in:
leonarski_f
2024-04-08 11:18:50 +02:00
parent 15d99c6162
commit c6d2b5eedf
72 changed files with 690 additions and 893 deletions
Download Patch File Download Diff File Expand all files Collapse all files

233
image_analysis/StrongPixelSet.cpp
View File

@@ -1,4 +1,10 @@
// Copyright (2019-2023) Paul Scherrer Institute
// Copyright (2019-2024) Paul Scherrer Institute
// SparseCCL code taken from https://github.com/acts-project/traccc/blob/main/core/include/traccc/clusterization/detail/sparse_ccl.hpp
// (c) 2021-2022 CERN for the benefit of the ACTS project
// Mozilla Public License Version 2.0
#include <bitset>
#include "StrongPixelSet.h"
@@ -57,145 +63,142 @@ void FilterSpotsByResolution(const DiffractionExperiment& experiment,
output.push_back(input[res_id_vector[i].id]);
}
StrongPixelSet::StrongPixelSet()
: strong_pixel_vector(RAW_MODULE_SIZE, false),
strong_pixel_count(0) {
StrongPixelSet::StrongPixelSet() : strong_pixel_count(0) {
pixels.reserve(max_strong_pixel_per_module);
}
void StrongPixelSet::AddStrongPixel(uint16_t col, uint16_t line, int32_t photons) {
auto key = strong_pixel_coord(col, line);
strong_pixel_map[key] = photons;
strong_pixel_vector.at(line * xpixel + col) = true;
pixels.push_back(strong_pixel{.col = col, .line = line, .counts = photons});
}
void StrongPixelSet::AddSingleStrongPixel(uint16_t col, uint16_t line, int32_t photons) {
auto key = strong_pixel_coord(col, line);
auto p = std::pair<uint32_t, int32_t>(key, photons);
single_pixels.emplace_back(p);
bool is_far_enough(strong_pixel pixel0, strong_pixel pixel1) {
return (pixel1.line - pixel0.line) > 1;
}
inline void StrongPixelSet::AddNeighbor(DiffractionSpot &spot, uint16_t col, uint16_t line) {
if (strong_pixel_vector[line * xpixel + col]) {
uint64_t coord = strong_pixel_coord(col, line);
auto iter = strong_pixel_map.find(coord);
ExtendSpot(spot, iter);
}
bool is_adjacent(strong_pixel pixel0, strong_pixel pixel1) {
return (fabs(pixel0.line - pixel1.line) <= 1) and
(fabs(pixel0.col - pixel1.col) <= 1);
}
// Creates a continuous spot
// strong pixels are loaded into dictionary (one dictionary per frame)
// and routine checks if neighboring pixels are also in dictionary (likely in log(N) time)
DiffractionSpot StrongPixelSet::BuildSpot(std::unordered_map<uint32_t, int32_t>::iterator &it) {
uint16_t col = col_from_strong_pixel(it->first);
uint16_t line = line_from_strong_pixel(it->first);
DiffractionSpot spot(col, line, it->second);
strong_pixel_vector[line * xpixel + col] = false;
strong_pixel_map.erase(it); // Remove strong pixel from the dictionary, so it is not processed again
if (col+1 < xpixel) {
AddNeighbor(spot, col + 1, line);
if (line < ypixel - 1)
AddNeighbor(spot, col + 1, line + 1);
if (line > 0)
AddNeighbor(spot, col + 1, line - 1);
uint16_t StrongPixelSet::find_root(uint16_t e) {
uint16_t r = e;
//assert(r < L.size());
while (L[r] != r) {
r = L[r];
//assert(r < L.size());
}
if (col != 0) {
AddNeighbor(spot, col - 1, line);
if (line < ypixel - 1)
AddNeighbor(spot, col - 1, line + 1);
if (line > 0)
AddNeighbor(spot, col - 1, line - 1);
}
if (line < ypixel - 1)
AddNeighbor(spot, col, line+1);
if (line > 0)
AddNeighbor(spot, col, line-1);
return spot;
return r;
}
void StrongPixelSet::ExtendSpot(DiffractionSpot &spot, std::unordered_map<uint32_t, int32_t>::iterator &it) {
uint16_t col = col_from_strong_pixel(it->first);
uint16_t line = line_from_strong_pixel(it->first);
spot.AddPixel(col, line, it->second);
strong_pixel_vector[line * xpixel + col] = false;
strong_pixel_map.erase(it); // Remove strong pixel from the dictionary, so it is not processed again
if (col+1 < xpixel) {
AddNeighbor(spot, col + 1, line);
if (line < ypixel - 1)
AddNeighbor(spot, col + 1, line + 1);
if (line > 0)
AddNeighbor(spot, col + 1, line - 1);
uint16_t StrongPixelSet::make_union(uint16_t e1, uint16_t e2) {
uint16_t e;
if (e1 < e2) {
e = e1;
//assert(e2 < L.size());
L[e2] = e;
} else {
e = e2;
//assert(e1 < L.size());
L[e1] = e;
}
if (col != 0) {
AddNeighbor(spot, col - 1, line);
if (line < ypixel - 1)
AddNeighbor(spot, col - 1, line + 1);
if (line > 0)
AddNeighbor(spot, col - 1, line - 1);
}
if (line < ypixel - 1)
AddNeighbor(spot, col, line+1);
if (line > 0)
AddNeighbor(spot, col, line-1);
return e;
}
void StrongPixelSet::FindSpots(const DiffractionExperiment &experiment, const SpotFindingSettings &settings,
std::vector<DiffractionSpot> &spots, uint16_t module_number) {
if (settings.min_pix_per_spot == 1) {
for (const auto &i: single_pixels) {
uint16_t col = col_from_strong_pixel(i.first);
uint16_t line = line_from_strong_pixel(i.first);
std::vector<DiffractionSpot> StrongPixelSet::sparseccl() {
L.resize(pixels.size());
DiffractionSpot spot(col, line, i.second);
spots.emplace_back(spot);
unsigned int labels = 0;
// first scan: pixel association
uint16_t start_j = 0;
for (uint16_t i = 0; i < pixels.size(); ++i) {
L[i] = i;
uint16_t ai = i;
for (uint16_t j = start_j; j < i; ++j) {
if (is_adjacent(pixels[i], pixels[j])) {
ai = make_union(ai, find_root(j));
} else if (is_far_enough(pixels[j], pixels[i])) {
++start_j;
}
}
}
while (!strong_pixel_map.empty()) {
auto iter = strong_pixel_map.begin();
DiffractionSpot spot = BuildSpot(iter);
spot.ConvertToImageCoordinates(experiment, module_number);
// second scan: transitive closure
for (unsigned int i = 0; i < L.size(); ++i) {
if (L[i] == i) {
L[i] = labels++;
} else {
L[i] = L[L[i]];
}
}
if ((spot.PixelCount() <= settings.max_pix_per_spot)
&& (spot.PixelCount() >= settings.min_pix_per_spot))
spots.push_back(spot);
std::vector<DiffractionSpot> spots(labels);
for (unsigned int i = 0; i < L.size(); i++)
spots[L[i]].AddPixel(pixels[i].col, pixels[i].line, pixels[i].counts);
return spots;
}
void StrongPixelSet::FindSpots(const DiffractionExperiment &experiment, const SpotFindingSettings &settings,
std::vector<DiffractionSpot> &spots, uint16_t module_number) {
if (!pixels.empty()) {
for (const auto &spot: sparseccl()) {
if ((spot.PixelCount() <= settings.max_pix_per_spot)
&& (spot.PixelCount() >= settings.min_pix_per_spot)) {
auto s = spot;
s.ConvertToImageCoordinates(experiment, module_number);
spots.push_back(s);
}
}
}
}
size_t StrongPixelSet::Count() const {
return strong_pixel_map.size();
}
size_t StrongPixelSet::Common(const StrongPixelSet &set) const {
size_t ret = 0;
for (const auto& pixel: strong_pixel_map) {
if (set.strong_pixel_map.count(pixel.first) == 1)
ret++;
}
return ret;
}
void StrongPixelSet::ReadFPGAOutput(const DeviceOutput &output) {
void StrongPixelSet::ReadFPGAOutput(const DiffractionExperiment & experiment,
const DeviceOutput &output) {
strong_pixel_count += output.spot_finding_result.strong_pixel_count;
for (int i = 0; i < std::min<uint32_t>(output.spot_finding_result.max_memory_index * 32,
SPOT_FINDER_MAX_STRONG_PIXEL); i++) {
uint32_t val = output.spot_finding_result.strong_pixel_number[i];
if (val != UINT32_MAX) {
uint32_t pixel = (val & ((1LU << 20) - 1));
uint8_t conn = ((val >> 24) & UINT8_MAX);
if (conn == 0)
AddSingleStrongPixel(pixel % RAW_MODULE_COLS, pixel / RAW_MODULE_COLS, output.pixels[pixel]);
else
AddStrongPixel(pixel % RAW_MODULE_COLS, pixel / RAW_MODULE_COLS, output.pixels[pixel]);
// Too many strong pixels will kill performance in data processing, so protection is needed
// Also if there are no strong pixels, there is no point in looking for them
if ((output.spot_finding_result.strong_pixel_count == 0) ||
(output.spot_finding_result.strong_pixel_count > max_strong_pixel_per_module))
return;
auto pixel_depth = experiment.GetPixelDepth();
auto out_ptr = (uint32_t *) output.spot_finding_result.strong_pixel;
if (pixel_depth == 2) {
for (int i = 0; i < RAW_MODULE_SIZE / (8 * sizeof(out_ptr[0])); i++) {
size_t npixel = i * 8 * sizeof(out_ptr[0]);
size_t line = npixel / RAW_MODULE_COLS;
if (out_ptr[i] != 0) {
std::bitset<32> bitset(out_ptr[i]);
for (int j = 0; j < 32; j++) {
if (bitset.test(j)) {
size_t col = (npixel | j) % RAW_MODULE_COLS;
AddStrongPixel(col, line, output.pixels[npixel]);
}
}
}
}
} else {
for (int i = 0; i < RAW_MODULE_SIZE / (8 * sizeof(out_ptr[0])); i++) {
size_t npixel = i * 8 * sizeof(out_ptr[0]);
size_t line = npixel / RAW_MODULE_COLS;
if (out_ptr[i] != 0) {
std::bitset<32> bitset(out_ptr[i]);
for (int j = 0; j < 32; j++) {
if (bitset.test(j)) {
size_t col = (npixel | j) % RAW_MODULE_COLS;
AddStrongPixel(col, line, ((int32_t *) output.pixels)[npixel]);
}
}
}
}
}
}