123 lines
4.7 KiB
C++
123 lines
4.7 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "BraggIntegrate2D.h"
|
|
#include "BraggPrediction.h"
|
|
|
|
#include "BraggIntegrationStats.h"
|
|
|
|
template<class T>
|
|
bool IntegrateReflection(Reflection &r, const T *image, size_t xpixel, size_t ypixel,
|
|
int64_t special_value, int64_t saturation,
|
|
float r_3, float r_1_sq, float r_2_sq, float r_3_sq) {
|
|
|
|
r.I = 0;
|
|
r.bkg = 0;
|
|
|
|
int64_t x0 = std::floor(r.center_x_pxl - r_3 - 1.0);
|
|
int64_t x1 = std::ceil(r.center_x_pxl + r_3 + 1.0);
|
|
int64_t y0 = std::floor(r.center_y_pxl - r_3 - 1.0);
|
|
int64_t y1 = std::ceil(r.center_y_pxl + r_3 + 1.0);
|
|
if (x0 < 0)
|
|
x0 = 0;
|
|
if (y0 < 0)
|
|
y0 = 0;
|
|
if (x1 >= xpixel)
|
|
x1 = xpixel - 1;
|
|
if (y1 >= ypixel)
|
|
y1 = ypixel - 1;
|
|
|
|
int64_t I_sum = 0;
|
|
int64_t I_npixel = 0;
|
|
int64_t bkg_sum = 0;
|
|
int64_t bkg_npixel = 0;
|
|
|
|
for (int64_t y = y0; y <= y1; y++) {
|
|
for (int64_t x = x0; x <= x1; x++) {
|
|
float dist_sq = (x - r.center_x_pxl) * (x - r.center_x_pxl) + (y - r.center_y_pxl) * (y - r.center_y_pxl);
|
|
if (image[y * xpixel + x] == special_value || image[y * xpixel + x] == saturation) {
|
|
continue;
|
|
} else if (dist_sq < r_1_sq) {
|
|
I_sum += image[y * xpixel + x];
|
|
I_npixel++;
|
|
} else if (dist_sq >= r_2_sq && dist_sq < r_3_sq) {
|
|
bkg_sum += image[y * xpixel + x];
|
|
bkg_npixel++;
|
|
}
|
|
}
|
|
}
|
|
|
|
if ((I_npixel > 2) && (bkg_npixel > 5)) {
|
|
r.bkg = static_cast<float>(bkg_sum) / static_cast<float>(bkg_npixel);
|
|
r.I = static_cast<float>(I_sum) - static_cast<float>(I_npixel) * r.bkg;
|
|
// minimum sigma is 1!
|
|
if (I_sum >= 1)
|
|
r.sigma = std::sqrt(static_cast<float>(I_sum));
|
|
else
|
|
r.sigma = 1;
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
|
|
template<class T>
|
|
QuickIntegrateResult IntegrateInternal(const DiffractionExperiment &experiment,
|
|
const CompressedImage &image,
|
|
const CrystalLattice &latt,
|
|
int64_t special_value,
|
|
int64_t saturation) {
|
|
auto settings = experiment.GetBraggIntegrationSettings();
|
|
|
|
|
|
auto start = std::chrono::high_resolution_clock::now();
|
|
std::vector<uint8_t> buffer;
|
|
auto ptr = reinterpret_cast<const T*>(image.GetUncompressedPtr(buffer));
|
|
|
|
auto refl = CalcBraggPredictions(experiment, latt, settings.GetDMinLimit_A(), settings.GetAngleFromEwald_deg());
|
|
|
|
QuickIntegrateResult ret;
|
|
BraggIntegrationStats stats(6.0, settings.GetDMinLimit_A(), 10); // Wilson statistics is only linear in certain resolution range
|
|
float r_3 = settings.GetR3();
|
|
float r_1_sq = settings.GetR1() * settings.GetR1();
|
|
float r_2_sq = settings.GetR2() * settings.GetR2();
|
|
float r_3_sq = settings.GetR3() * settings.GetR3();
|
|
|
|
for (const auto &[x, r_in]: refl) {
|
|
Reflection r = r_in;
|
|
if (IntegrateReflection(r, ptr, image.GetWidth(), image.GetHeight(), special_value, saturation,
|
|
r_3, r_1_sq, r_2_sq, r_3_sq)) {
|
|
ret.reflections.push_back(r);
|
|
stats.AddReflection(r);
|
|
}
|
|
}
|
|
ret.b_factor = stats.BFactor();
|
|
return ret;
|
|
}
|
|
|
|
QuickIntegrateResult BraggIntegrate2D(const DiffractionExperiment &experiment,
|
|
const CompressedImage &image,
|
|
const CrystalLattice &latt) {
|
|
if (image.GetCompressedSize() == 0)
|
|
return {};
|
|
|
|
switch (image.GetMode()) {
|
|
case CompressedImageMode::Int8:
|
|
return IntegrateInternal<int8_t>(experiment, image, latt, INT8_MIN, INT8_MAX);
|
|
case CompressedImageMode::Int16:
|
|
return IntegrateInternal<int16_t>(experiment, image, latt, INT16_MIN, INT16_MAX);
|
|
case CompressedImageMode::Int32:
|
|
return IntegrateInternal<int32_t>(experiment, image, latt, INT32_MIN, INT32_MAX);
|
|
case CompressedImageMode::Uint8:
|
|
return IntegrateInternal<uint8_t>(experiment, image, latt, UINT8_MAX, UINT8_MAX);
|
|
case CompressedImageMode::Uint16:
|
|
return IntegrateInternal<uint16_t>(experiment, image, latt, UINT16_MAX, UINT16_MAX);
|
|
case CompressedImageMode::Uint32:
|
|
return IntegrateInternal<uint16_t>(experiment, image, latt, UINT32_MAX, UINT32_MAX);
|
|
default:
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Image mode not supported");
|
|
}
|
|
}
|
|
|