Commit ba4f4fe7 authored by Alejandro Homs Puron's avatar Alejandro Homs Puron Committed by Generic Bliss account for Control Software
Browse files

[PROC][FAI] Add average of corrected and peak found images

parent 7cf850fc
......@@ -73,6 +73,7 @@ namespace processing::fai
vector<float> radius1d; //
float noise; //
float cutoff_pick; //
int acc_update_freq; // Accumulation
};
namespace detail
......@@ -129,6 +130,7 @@ namespace processing::fai
std::filesystem::path radius2d_path, // Peak finding
std::filesystem::path radius1d_path, //
float noise, float cutoff_pick, //
int acc_update_freq, // Accumulation
boost::compute::context& context, boost::compute::command_queue& queue)
{
kernel_params res(context);
......@@ -154,6 +156,7 @@ namespace processing::fai
res.empty = empty;
res.noise = noise;
res.cutoff_pick = cutoff_pick;
res.acc_update_freq = acc_update_freq;
if (!mask_path.empty())
detail::read_h5_dset(mask_path, res.mask, "/data", queue);
......@@ -231,6 +234,8 @@ namespace processing::fai
extract(noise, "noise");
float cutoff_pick;
extract(cutoff_pick, "cutoff_pick");
int acc_update_freq;
extract(acc_update_freq, "acc_update_freq");
// Preprocessing
std::filesystem::path variance_path;
std::string mask_path;
......@@ -277,6 +282,7 @@ namespace processing::fai
radius2d_path, // Peak finding
radius1d_path, //
noise, cutoff_pick, //
acc_update_freq, //
context, queue);
}
......
......@@ -72,6 +72,7 @@ namespace processing::fai
vector<float> const& radius1d, //
float noise, //
float cutoff_pick, //
int acc_update_freq, // Accumulation
std::filesystem::path const& cl_source_path, // Path to the CL sources
bcl::command_queue& queue) //
{
......@@ -115,6 +116,7 @@ namespace processing::fai
bool do_absorption = !absorption_d.empty();
bool do_mask = !mask_d.empty();
bool do_dummy = false;
bool do_acc = (acc_update_freq > 0);
bcl::vector<bcl::float4_> preprocessed_d(nb_pixels, context); //!< Preprocessed data
bcl::copy(mask.begin(), mask.end(), mask_d.begin(), queue);
bcl::kernel stage2(stage2_program, "corrections4");
......@@ -141,39 +143,51 @@ namespace processing::fai
bcl::vector<int> counter_d(1, context);
bcl::kernel stage4(stage4_program, "find_peaks");
// Stage 6 - Accumulate
auto acc_nb_pixels = acc_update_freq ? nb_pixels : 0;
bcl::vector<float> acc_corr_d(acc_nb_pixels, 0.0f, queue); //!< Accum pixels after correction
bcl::vector<float> acc_peak_d(acc_nb_pixels, 0.0f, queue); //!< Accum pixels after peak find
std::uint32_t curr_acc_frames = 0;
const bool reset_acc = false;
return [=, //
gain = std::move(gain_d),
pedestal = std::move(pedestal_d), //
corrected = std::move(corrected_d), //
//
variance = std::move(variance_d), //
dark = std::move(dark_d), //
dark_variance = std::move(dark_variance_d), //
flat = std::move(flat_d), //
solid_angle = std::move(solid_angle_d), //
polarization = std::move(polarization), //
absorption = std::move(absorption_d), //
mask = std::move(mask_d), //
preprocessed = std::move(preprocessed_d), //
//
csr_data = std::move(csr_data_d), //
csr_indices = std::move(csr_indices_d), //
csr_indptr = std::move(csr_indptr_d), //
summed = std::move(summed_d), //
stderrmean = std::move(stderrmean_d), //
//
radius2d = std::move(radius2d_d), //
radius1d = std::move(radius1d_d), //
counter = std::move(counter_d) //
](bcl::command_queue& queue, //
vector<std::uint16_t> const& raw, // Input
vector<float>& corrected_dbg, // Corrected
vector<bcl::float4_>& preprocessed_dbg, // Preprocessed
vector<bcl::float8_>& cliped_dbg, // Cliped (after sigma clip)
vector<bcl::float4_>& found_dbg, // Found (after peak find)
vector<int>& peak_idx, vector<float>& peak_val, // Outputs Peaks
vector<float>& bkg_avg, vector<float>& bkg_std // Background
) mutable {
variance = std::move(variance_d), //
dark = std::move(dark_d), //
dark_variance = std::move(dark_variance_d), //
flat = std::move(flat_d), //
solid_angle = std::move(solid_angle_d), //
polarization = std::move(polarization), //
absorption = std::move(absorption_d), //
mask = std::move(mask_d), //
preprocessed = std::move(preprocessed_d), //
//
csr_data = std::move(csr_data_d), //
csr_indices = std::move(csr_indices_d), //
csr_indptr = std::move(csr_indptr_d), //
summed = std::move(summed_d), //
stderrmean = std::move(stderrmean_d), //
//
radius2d = std::move(radius2d_d), //
radius1d = std::move(radius1d_d), //
counter = std::move(counter_d), //
//
acc_corr_d = std::move(acc_corr_d), //
acc_peak_d = std::move(acc_peak_d) //
](bcl::command_queue& queue, //
vector<std::uint16_t> const& raw, // Input
vector<float>& corrected_dbg, // Corrected
vector<bcl::float4_>& preprocessed_dbg, // Preprocessed
vector<bcl::float8_>& cliped_dbg, // Cliped (after sigma clip)
vector<bcl::float4_>& found_dbg, // Found (after peak find)
vector<int>& peak_idx, vector<float>& peak_val, // Outputs Peaks
vector<float>& bkg_avg, vector<float>& bkg_std, // Background
vector<float>& acc_corr, vector<float>& acc_peak, // Accumulation
bool force_acc) mutable {
// Stage 0 - Initialization (might not be necessary)
stage0.set_args(bkg_avg, bkg_std, summed, (unsigned int) nb_bins);
queue.enqueue_1d_range_kernel(stage0, 0, nb_bins, 0 /*max = 1024 for Tesla*/);
......@@ -191,7 +205,7 @@ namespace processing::fai
(char) do_dark_variance, dark_variance, (char) do_flat, flat, (char) do_solid_angle,
solid_angle, (char) do_polarization, polarization, (char) do_absorption, absorption,
(char) do_mask, mask, (char) do_dummy, dummy, delta_dummy, normalization_factor,
preprocessed, (unsigned int) nb_pixels);
(char) do_acc, acc_corr_d, preprocessed, (unsigned int) nb_pixels);
queue.enqueue_1d_range_kernel(stage2, 0, nb_pixels, 0);
// DBG
......@@ -213,7 +227,8 @@ namespace processing::fai
// Stage 4 - Peak finding
bcl::fill(counter.begin(), counter.end(), 0, queue);
stage4.set_args(preprocessed, radius2d, radius1d, bkg_avg, bkg_std, radius_min, radius_max, cutoff_pick,
noise, (unsigned int) nb_pixels, (unsigned int) nb_bins, counter, peak_idx);
noise, (unsigned int) nb_pixels, (unsigned int) nb_bins, counter, peak_idx, (char) do_acc,
acc_peak_d);
auto nb_groups = (nb_pixels - 1) / peak_finder_wg_size + 1;
queue.enqueue_1d_range_kernel(stage4, 0, nb_groups * peak_finder_wg_size, peak_finder_wg_size);
......@@ -229,20 +244,38 @@ namespace processing::fai
{ return preprocessed[idx].s1; });
bcl::transform(peak_idx.begin(), peak_idx.begin() + nb_peaks, peak_val.begin(), copy_peaks, queue);
return nb_peaks;
// Stage 6 - Accumulate corr & peaks
std::size_t ret_acc_frames = 0;
if (do_acc) {
// Increment nb of accumulated samples and check if update is needed
if ((++curr_acc_frames % acc_update_freq == 0) || force_acc) {
ret_acc_frames = curr_acc_frames;
bcl::copy(acc_corr_d.begin(), acc_corr_d.end(), acc_corr.begin(), queue);
bcl::copy(acc_peak_d.begin(), acc_peak_d.end(), acc_peak.begin(), queue);
if (reset_acc) {
curr_acc_frames = 0;
bcl::fill(acc_corr_d.begin(), acc_corr_d.end(), 0.0f, queue);
bcl::fill(acc_peak_d.begin(), acc_peak_d.end(), 0.0f, queue);
}
}
}
return std::make_tuple(std::size_t(nb_peaks), ret_acc_frames);
};
};
//using peak_finder = decltype(make_peak_finder(std::declval<bcl::context>(), std::declval<bcl::command_queue>()));
using peak_finder = std::function<int(bcl::command_queue&, //
vector<std::uint16_t> const&, // Input
vector<float>& corrected_dbg, // Corrected
vector<bcl::float4_>& preprocessed_dbg, // Preprocessed
vector<bcl::float8_>& cliped_dbg, // Summed (after sigma clip)
vector<bcl::float4_>& found_dbg, // Finded (after peak find)
vector<int>&, vector<float>&, // Outputs Peaks
vector<float>&, vector<float>& // Background
)>;
using ret_t = std::tuple<std::size_t, std::size_t>;
using peak_finder = std::function<ret_t(bcl::command_queue&, //
vector<std::uint16_t> const&, // Input
vector<float>& corrected_dbg, // Corrected
vector<bcl::float4_>& preprocessed_dbg, // Preprocessed
vector<bcl::float8_>& cliped_dbg, // Summed (after sigma clip)
vector<bcl::float4_>& found_dbg, // Finded (after peak find)
vector<int>&, vector<float>&, // Outputs Peaks
vector<float>&, vector<float>&, // Background
vector<float>&, vector<float>&, bool // Accumulation
)>;
} // namespace processing::fai
} // namespace lima
......@@ -23,7 +23,9 @@ __kernel void find_peaks( global float4 *preproc4, // both input and
const unsigned int nb_pixels,
const unsigned int nb_bins,
global int *counter, // Counter of the number of peaks found
global int *highidx // indexes of the pixels of high intensity
global int *highidx, // indexes of the pixels of high intensity
const char do_acc, // Perform accumulation
global float *acc_peak // Accumulated peaks
){
int tid = get_local_id(0);
int gid = get_global_id(0);
......@@ -66,6 +68,8 @@ __kernel void find_peaks( global float4 *preproc4, // both input and
value.s3 = fast_length((float2)(value.s3, noise)); //add quadratically noise to std
if ((value.s1 - value.s2) > max(noise, cutoff*value.s3)){
local_highidx[atomic_inc(local_counter)] = gid;
if (do_acc)
acc_peak[gid] += value.s1;
}//pixel is considered of high intensity: registering it.
} //check radius range
preproc4[gid] = value;
......
......@@ -32,7 +32,8 @@ __kernel void corrections4(const __global float* in, const char poissonian, cons
const __global float* polarization, const char do_absorption,
const __global float* absorption, const char do_mask, const __global char* mask,
const char do_dummy, const float dummy, const float delta_dummy,
const float normalization_factor, __global float4* res, const unsigned int length)
const float normalization_factor, const char do_acc, __global float* acc,
__global float4* res, const unsigned int length)
{
size_t i = get_global_id(0);
float4 result = (float4)(in[i], 0.0f, 0.0f, 0.0f);
......@@ -65,6 +66,8 @@ __kernel void corrections4(const __global float* in, const char poissonian, cons
result.s2 *= absorption[i];
if (isnan(result.s0) || isnan(result.s1) || isnan(result.s2) || (result.s2 == 0.0f))
result = (float4)(0.0f, 0.0f, 0.0f, 0.0f);
else if (do_acc)
acc[i] += result.s0;
} else {
result = (float4)(0.0f, 0.0f, 0.0f, 0.0f);
}
......
......@@ -111,6 +111,7 @@ struct peak_finder
py::c_array_t<float> radius1d_np, //
float noise, //
float cutoff_pick, //
int acc_update_freq, //
std::string const& cl_source_path // CL source file path
) :
context(context_arg),
......@@ -122,7 +123,9 @@ struct peak_finder
peak_indices_d(width * height, context_arg),
peak_values_d(width * height, context_arg),
background_avg_d(nb_bins, context_arg),
background_std_d(nb_bins, context_arg)
background_std_d(nb_bins, context_arg),
acc_corr_d(width * height, context_arg),
acc_peak_d(width * height, context_arg)
{
// Create queue
queue = bcl::command_queue{context, context.get_device()};
......@@ -178,20 +181,22 @@ struct peak_finder
radius2d, radius1d, // Peak finding
noise, //
cutoff_pick, //
acc_update_freq, //
cl_source_path, // CL source file path
queue);
}
~peak_finder() { queue.finish(); }
auto apply(py::c_array_t<std::uint16_t> raw)
auto apply(py::c_array_t<std::uint16_t> raw, bool force_acc = false)
{
// Write input to the device (non-blocking)
copy_async(raw, raw_d, queue);
// Run the kernels
std::size_t nb_peaks = pf(queue, raw_d, corrected_dbg_d, preprocessed_dbg_d, cliped_dbg_d, found_dbg_d,
peak_indices_d, peak_values_d, background_avg_d, background_std_d);
auto [nb_peaks, nb_acc_frames] =
pf(queue, raw_d, corrected_dbg_d, preprocessed_dbg_d, cliped_dbg_d, found_dbg_d, peak_indices_d,
peak_values_d, background_avg_d, background_std_d, acc_corr_d, acc_peak_d, force_acc);
py::c_array_t<float> corrected({(int) corrected_dbg_d.size()});
py::c_array_t<float> preprocessed({(int) preprocessed_dbg_d.size(), 4});
......@@ -202,12 +207,18 @@ struct peak_finder
py::c_array_t<float> peak_values(nb_peaks);
py::c_array_t<float> background_avg(background_avg_d.size());
py::c_array_t<float> background_std(background_std_d.size());
py::c_array_t<float> acc_corr(nb_acc_frames ? acc_corr_d.size() : 0);
py::c_array_t<float> acc_peak(nb_acc_frames ? acc_peak_d.size() : 0);
// Read result from the device to array (non-blocking)
copy_n(peak_indices_d, nb_peaks, peak_indices, queue);
copy_n(peak_values_d, nb_peaks, peak_values, queue);
copy_async(background_avg_d, background_avg, queue);
copy_async(background_std_d, background_std, queue);
if (nb_acc_frames) {
copy_async(acc_corr_d, acc_corr, queue);
copy_async(acc_peak_d, acc_peak, queue);
}
copy_async(corrected_dbg_d, corrected, queue);
copy_async(preprocessed_dbg_d, preprocessed, queue);
......@@ -215,7 +226,7 @@ struct peak_finder
copy_async(found_dbg_d, found, queue);
return std::make_tuple(corrected, preprocessed, cliped, found, peak_indices, peak_values, background_avg,
background_std);
background_std, nb_acc_frames, acc_corr, acc_peak);
}
fai::vector<std::uint16_t> raw_d;
......@@ -227,6 +238,8 @@ struct peak_finder
fai::vector<float> peak_values_d;
fai::vector<float> background_avg_d;
fai::vector<float> background_std_d;
fai::vector<float> acc_corr_d;
fai::vector<float> acc_peak_d;
bcl::context context;
bcl::command_queue queue;
......@@ -288,6 +301,7 @@ PYBIND11_MODULE(peak_finder, m)
py::c_array_t<float> radius1d, //
float noise, //
float cutoff_pick, //
int acc_update_freq, //
std::string const& cl_source_path // Path to the CL source files
) {
// Filter for a 1.2 platform and set it as the default
......@@ -350,6 +364,7 @@ PYBIND11_MODULE(peak_finder, m)
radius2d, radius1d, // Peak finding
noise, //
cutoff_pick, //
acc_update_freq, //
cl_source_path); // CL source file path
}),
"Construct a PeakFinder callable", //
......@@ -358,9 +373,10 @@ PYBIND11_MODULE(peak_finder, m)
"dark"_a, "dark_variance"_a, "flat"_a, "solidangle"_a, "polarization"_a, "absorption"_a,
"normalization_factor"_a = 1.0, "csr_data"_a, "csr_indices"_a, "csr_indptr"_a, "cutoff_clip"_a = 5.0,
"cycle"_a = 5, "empty"_a = -1.0, "radius2d"_a, "radius1d"_a, "noise"_a = 0.5, "cutoff_pick"_a = 3.0,
"cl_source_path"_a = "detectors/psi/processing/src")
"acc_update_freq"_a = 0, "cl_source_path"_a = "detectors/psi/processing/src")
.def("__call__", &pyfai::peak_finder::apply,
"Perform a sigma-clipping iterative filter within each bin along each row.", "data"_a)
"Perform a sigma-clipping iterative filter within each bin along each row.", "data"_a,
"force_acc"_a = false)
.def("__repr__", [](const pyfai::peak_finder& pf) {
return fmt::format("<PeakFinder on '{}'>", pf.context.get_device().name());
});
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment