Jungfraujoch/image_analysis/ImageSpotFinder.cpp

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

#include "ImageSpotFinder.h"
#include "../common/JFJochException.h"

ImageSpotFinder::ImageSpotFinder(const DiffractionGeometry &geom, size_t width, size_t height)
: resolution_map(width*height), width(width), height(height) {
    for (int row = 0; row < height; row++) {
        for (int col = 0; col < width; col++) {
            resolution_map[row * width + col] = geom.PxlToRes(col, row);
        }
    }
}

ImageSpotFinder::ImageSpotFinder(std::vector<float> &in_resolution_map, size_t in_width, size_t in_height)
: resolution_map(in_resolution_map), width(in_width), height(in_height) {
    if (resolution_map.size() != width * height)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Coordinate mismatch");
}

ImageSpotFinder &ImageSpotFinder::PixelMask(std::vector<uint32_t> &pixel_mask) {
    if (pixel_mask.size() != width * height)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "Coordinate mismatch");

    for (int pxl = 0; pxl < width *height; pxl++) {
        if (pixel_mask[pxl] != 0)
            resolution_map[pxl] = 0.0f;
    }
    return *this;
}

void ImageSpotFinder::nbx(int input) {
    NBX = input;
}

std::vector<DiffractionSpot> ImageSpotFinder::Run(const int32_t *input,
                                                  const SpotFindingSettings &settings) const {
    std::vector<DiffractionSpot> output;
    StrongPixelSet strong_pixel_set;

    if (settings.signal_to_noise_threshold <= 0.0) {
        if (settings.photon_count_threshold > 0) {
            for (int line = 0; line < height; line++) {
                for (int col = 0; col < width; col++) {
                    int32_t pxl_val = input[line * width + col];
                    if (pxl_val > settings.photon_count_threshold
                        && (resolution_map[line * width + col] >= settings.high_resolution_limit)
                        && (resolution_map[line * width + col] <= settings.low_resolution_limit))
                        strong_pixel_set.AddStrongPixel(col, line, pxl_val);
                }
            }
            strong_pixel_set.FindSpotsImage(settings, output);
        }
        return output;
    }

    float strong2 = settings.signal_to_noise_threshold * settings.signal_to_noise_threshold;

    std::vector<uint8_t> mask(width * height, 0);

    // Sum and sum of squares of (2*NBY+1) vertical elements
    // These are updated after each line is finished
    // 64-bit integer guarantees calculations are made without rounding errors
    std::vector<int64_t> sum_vert(width, 0);
    std::vector<int64_t> sum2_vert(width, 0);
    std::vector<uint16_t> valid_vert(width, 0);

    for (int line = 0; line < NBX; line++) {
        for (int col = 0; col < width; col++) {
            auto pxl = line * width + col;
            mask[pxl] = (input[pxl] == INT32_MIN)
                        || (resolution_map[pxl] < settings.high_resolution_limit)
                        || (resolution_map[pxl] > settings.low_resolution_limit);

            if (mask[pxl] == 0) {
                int64_t tmp = input[pxl];
                sum_vert[col] += tmp;
                sum2_vert[col] += tmp * tmp;
                valid_vert[col] += 1;
            }
        }
    }

    for (int line = 0; line < height; line++) {
        for (int col = 0; col < width; col++) {
            if (line < height - NBX) {
                auto pxl = (line + NBX) * width + col;
                mask[pxl] = (input[pxl] == INT32_MIN)
                            || (resolution_map[pxl] <= settings.high_resolution_limit)
                            || (resolution_map[pxl] >= settings.low_resolution_limit);

                if (mask[pxl] == 0) {
                    int64_t tmp = input[pxl];
                    sum_vert[col] += tmp;
                    sum2_vert[col] += tmp * tmp;
                    valid_vert[col] += 1;
                }
            }

            if (line >= NBX + 1) {
                auto pxl = (line - (NBX + 1)) * width + col;
                if (mask[pxl] == 0) {
                    int64_t tmp = input[pxl];
                    sum_vert[col] -= tmp;
                    sum2_vert[col] -= tmp * tmp;
                    valid_vert[col] -= 1;
                }
            }
        }

        int64_t sum = 0;
        int64_t sum2 = 0;
        int64_t valid = 0;

        for (int col = 0; col < NBX; col++) {
            sum += sum_vert[col];
            sum2 += sum2_vert[col];
            valid += valid_vert[col];
        }

        for (int col = 0; col < width; col++) {
            if (col < width - NBX) {
                sum += sum_vert[col + NBX];
                sum2 += sum2_vert[col + NBX];
                valid += valid_vert[col + NBX];
            }

            if (col >= NBX + 1) {
                sum -= sum_vert[col - NBX - 1];
                sum2 -= sum2_vert[col - NBX - 1];
                valid -= valid_vert[col - NBX - 1];
            }

            int32_t pxl_val = input[line * width + col];
            int64_t var = valid * sum2 - (sum * sum);
            int64_t in_minus_mean = pxl_val * valid - sum;

            if ((pxl_val > settings.photon_count_threshold) && // pixel is above count threshold
                (in_minus_mean > 0) && // pixel value is larger than mean
                (in_minus_mean * in_minus_mean > std::ceil(var * strong2)) && // pixel is above SNR threshold
                (mask[line * width + col] == 0)) // pixel is not bad pixel
                strong_pixel_set.AddStrongPixel(col, line, pxl_val);
        }
    }

    strong_pixel_set.FindSpotsImage(settings, output);
    return output;
}