Jungfraujoch/image_analysis/spot_finding/StrongPixelSet.cpp

// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

// 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"

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) {
    pixels.push_back(strong_pixel{.col = col, .line = line, .counts = photons});
    ++strong_pixel_count;
}

bool is_far_enough(strong_pixel pixel0, strong_pixel pixel1) {
    return (pixel1.line - pixel0.line) > 1;
}

bool is_adjacent(strong_pixel pixel0, strong_pixel pixel1) {
    auto line_diff = pixel0.line - pixel1.line;
    auto col_diff = pixel0.col - pixel1.col;
    return line_diff <= 1 && line_diff >= -1 && col_diff <= 1 && col_diff >= -1;
}

uint32_t StrongPixelSet::find_root(uint32_t e) {
    uint32_t r = e;
    while (L[r] != r)
        r = L[r];
    return r;
}

uint32_t StrongPixelSet::make_union(uint32_t e1, uint32_t e2) {
    uint32_t e;
    if (e1 < e2) {
        e = e1;
        L[e2] = e;
    } else {
        e = e2;
        L[e1] = e;
    }
    return e;
}

std::vector<DiffractionSpot> StrongPixelSet::sparseccl() {
    L.resize(pixels.size());

    unsigned int labels = 0;

    // first scan: pixel association
    uint32_t start_j = 0;
    for (uint32_t i = 0; i < pixels.size(); ++i) {
        L[i] = i;
        uint32_t ai = i;
        for (uint32_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;
            }
        }
    }

    // second scan: transitive closure
    for (uint32_t i = 0; i < L.size(); ++i) {
        if (L[i] == i) {
            L[i] = labels++;
        } else {
            L[i] = L[L[i]];
        }
    }

    std::vector<DiffractionSpot> spots(labels);

    for (uint32_t i = 0; i < L.size(); i++)
        spots[L[i]].AddPixel(pixels[i].col, pixels[i].line, pixels[i].counts);

    return spots;
}


void StrongPixelSet::FindSpotsImage(const SpotFindingSettings &settings, std::vector<DiffractionSpot> &spots) {
    // Avoid spot finding, when more than 65536 strong pixel count (will be super slow)
    if (!pixels.empty() && (strong_pixel_count < UINT16_MAX)) {
        for (const auto &spot: sparseccl()) {
            if ((spot.PixelCount() <= settings.max_pix_per_spot)
                && (spot.PixelCount() >= settings.min_pix_per_spot)) {
                spots.push_back(spot);
            }
        }
    }
}

void StrongPixelSet::FindSpots(const DiffractionExperiment &experiment, const SpotFindingSettings &settings,
                               std::vector<DiffractionSpot> &spots, uint16_t module_number) {
    // Avoid spot finding, when more than 65536 strong pixel count (will be super slow)
    if (!pixels.empty() && (strong_pixel_count < UINT16_MAX)) {
        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);
            }
        }
    }
}

void StrongPixelSet::ReadFPGAOutput(const DiffractionExperiment & experiment,
                                    const DeviceOutput &output) {
    // 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)) {
        // If max strong pixel per module condition kicks-in, still report correct strong pixel count
        strong_pixel_count = output.spot_finding_result.strong_pixel_count;
        return;
    }

    auto pixel_depth = experiment.GetByteDepthImage();
    auto out_ptr = (uint32_t *) output.spot_finding_result.strong_pixel;

    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;
                    if (pixel_depth == 2)
                        AddStrongPixel(col, line, output.pixels[npixel | j]);
                    else if (pixel_depth == 1)
                        AddStrongPixel(col, line, ((int8_t *)output.pixels)[npixel | j]);
                    else if (pixel_depth == 4)
                        AddStrongPixel(col, line, ((int32_t *)output.pixels)[npixel | j]);
                }
            }
        }
    }

}

uint32_t StrongPixelSet::GetStrongPixelCount() const {
    return strong_pixel_count;
}