diff --git a/setup.py b/setup.py index 21caef8e7f104d4eb63837f0d14f1590c684774e..c9555f58b91f7797dbefee04fff01c59caca9bbb 100644 --- a/setup.py +++ b/setup.py @@ -30,6 +30,7 @@ setup(name='calng', 'ManualDsscGeometry = calng.DsscCorrection:ManualDsscGeometry', 'JungfrauCorrection = calng.JungfrauCorrection:JungfrauCorrection', 'ManualJungfrauGeometry = calng.JungfrauCorrection:ManualJungfrauGeometry', + 'LpdCorrection = calng.LpdCorrection:LpdCorrection', 'ShmemToZMQ = calng.ShmemToZMQ:ShmemToZMQ', 'ShmemTrainMatcher = calng.ShmemTrainMatcher:ShmemTrainMatcher', 'DetectorAssembler = calng.DetectorAssembler:DetectorAssembler', @@ -38,6 +39,8 @@ setup(name='calng', 'karabo.middlelayer_device': [ 'CalibrationManager = calng.CalibrationManager:CalibrationManager', 'MdlAgipd1MGeometry = calng.mdl_geometry_base:MdlAgipd1MGeometry', + 'MdlDssc1MGeometry = calng.mdl_geometry_base:MdlDssc1MGeometry', + 'MdlLpd1MGeometry = calng.mdl_geometry_base:MdlLpd1MGeometry', 'MdlJungfrauGeometry = calng.mdl_geometry_base:MdlJungfrauGeometry', ], }, diff --git a/src/calng/AgipdCorrection.py b/src/calng/AgipdCorrection.py index 353de48cb44ba2d75325cdf7b242ad23582870fc..caa7812efe602e8f0457afdbbdedd933262e24fa 100644 --- a/src/calng/AgipdCorrection.py +++ b/src/calng/AgipdCorrection.py @@ -36,35 +36,6 @@ class AgipdGainMode(enum.IntEnum): FIXED_MEDIUM_GAIN = 2 FIXED_LOW_GAIN = 3 - -class BadPixelValues(enum.IntFlag): - """The European XFEL Bad Pixel Encoding - - Straight from pycalibration's enum.py""" - - OFFSET_OUT_OF_THRESHOLD = 2 ** 0 - NOISE_OUT_OF_THRESHOLD = 2 ** 1 - OFFSET_NOISE_EVAL_ERROR = 2 ** 2 - NO_DARK_DATA = 2 ** 3 - CI_GAIN_OUT_OF_THRESHOLD = 2 ** 4 - CI_LINEAR_DEVIATION = 2 ** 5 - CI_EVAL_ERROR = 2 ** 6 - FF_GAIN_EVAL_ERROR = 2 ** 7 - FF_GAIN_DEVIATION = 2 ** 8 - FF_NO_ENTRIES = 2 ** 9 - CI2_EVAL_ERROR = 2 ** 10 - VALUE_IS_NAN = 2 ** 11 - VALUE_OUT_OF_RANGE = 2 ** 12 - GAIN_THRESHOLDING_ERROR = 2 ** 13 - DATA_STD_IS_ZERO = 2 ** 14 - ASIC_STD_BELOW_NOISE = 2 ** 15 - INTERPOLATED = 2 ** 16 - NOISY_ADC = 2 ** 17 - OVERSCAN = 2 ** 18 - NON_SENSITIVE = 2 ** 19 - NON_LIN_RESPONSE_REGION = 2 ** 20 - - class CorrectionFlags(enum.IntFlag): NONE = 0 THRESHOLD = 1 @@ -564,7 +535,7 @@ class AgipdCorrection(BaseCorrection): AgipdCorrection._managed_keys.add("corrections.gainXray.gGainValue") AgipdCorrection._managed_keys.add("corrections.badPixels.maskingValue") # TODO: DRY / encapsulate - for field in BadPixelValues: + for field in utils.BadPixelValues: ( BOOL_ELEMENT(expected) .key(f"corrections.badPixels.subsetToUse.{field.name}") @@ -752,7 +723,7 @@ class AgipdCorrection(BaseCorrection): def _update_bad_pixel_selection(self): selection = 0 - for field in BadPixelValues: + for field in utils.BadPixelValues: if self.get(f"corrections.badPixels.subsetToUse.{field.name}"): selection |= field self._override_bad_pixel_flags = selection @@ -816,7 +787,7 @@ class AgipdCorrection(BaseCorrection): update.get( f"corrections.badPixels.subsetToUse.{field.name}", default=False ) - for field in BadPixelValues + for field in utils.BadPixelValues ): self.log_status_info( "Some fields reenabled, reloading cached bad pixel constants" diff --git a/src/calng/DsscCorrection.py b/src/calng/DsscCorrection.py index 865660755d9179318b649d0527c294ff4aa04f17..7bf97300b86e63f6295127e15dd8657df607724b 100644 --- a/src/calng/DsscCorrection.py +++ b/src/calng/DsscCorrection.py @@ -219,7 +219,7 @@ class DsscCorrection(BaseCorrection): cell_table, do_generate_preview, ): - pulse_table = np.squeeze(data_hash.get("image.pulseId")) + pulse_table = np.ravel(data_hash.get("image.pulseId")) if self._frame_filter is not None: try: cell_table = cell_table[self._frame_filter] diff --git a/src/calng/LpdCorrection.py b/src/calng/LpdCorrection.py new file mode 100644 index 0000000000000000000000000000000000000000..6d2505d8adf6a58f1ee5cb4ad1946276eb1dba77 --- /dev/null +++ b/src/calng/LpdCorrection.py @@ -0,0 +1,401 @@ +import enum + +import cupy +import numpy as np +from karabo.bound import ( + DOUBLE_ELEMENT, + KARABO_CLASSINFO, + OVERWRITE_ELEMENT, + VECTOR_STRING_ELEMENT, +) + +from . import base_gpu, base_calcat, utils +from ._version import version as deviceVersion +from .base_correction import BaseCorrection, add_correction_step_schema + + +class CorrectionFlags(enum.IntFlag): + NONE = 0 + OFFSET = 1 + GAIN_AMP = 2 + REL_GAIN = 4 + FF_CORR = 8 + BPMASK = 16 + + +class LpdConstants(enum.Enum): + Offset = enum.auto() + BadPixelsDark = enum.auto() + GainAmpMap = enum.auto() + FFMap = enum.auto() + RelativeGain = enum.auto() + BadPixelsFF = enum.auto() + + +class LpdGpuRunner(base_gpu.BaseGpuRunner): + _kernel_source_filename = "lpd_gpu.cu" + _corrected_axis_order = "cxy" # TODO: get specs for LPD data + + def __init__( + self, + pixels_x, + pixels_y, + memory_cells, + constant_memory_cells, + input_data_dtype=cupy.uint16, + output_data_dtype=cupy.float32, + bad_pixel_mask_value=cupy.nan, + ): + self.input_shape = (memory_cells, 1, pixels_y, pixels_x) + self.processed_shape = (memory_cells, pixels_y, pixels_x) + super().__init__( + pixels_x, + pixels_y, + memory_cells, + constant_memory_cells, + input_data_dtype, + output_data_dtype, + ) + self.map_shape = (constant_memory_cells, pixels_y, pixels_x, 3) + self.offset_map_gpu = cupy.zeros(self.map_shape, dtype=cupy.float32) + self.gain_amp_map_gpu = cupy.ones(self.map_shape, dtype=cupy.float32) + self.rel_gain_intercept_map_gpu = cupy.zeros(self.map_shape, dtype=cupy.float32) + self.rel_gain_slopes_map_gpu = cupy.ones(self.map_shape, dtype=cupy.float32) + self.flatfield_map_gpu = cupy.ones(self.map_shape, dtype=cupy.float32) + self.bad_pixel_map_gpu = cupy.ones(self.map_shape, dtype=cupy.uint32) + self.bad_pixel_mask_value = bad_pixel_mask_value + + self.update_block_size((1, 1, 64)) + + def _get_raw_for_preview(self): + return self.input_data_gpu[:, 0] + + def _get_corrected_for_preview(self): + return self.processed_data_gpu + + def correct(self, flags): + self.correction_kernel( + self.full_grid, + self.full_block, + ( + self.input_data_gpu, + self.cell_table_gpu, + cupy.uint8(flags), + self.offset_map_gpu, + self.gain_amp_map_gpu, + self.rel_gain_intercept_map_gpu, + self.rel_gain_slopes_map_gpu, + self.flatfield_map_gpu, + self.bad_pixel_map_gpu, + self.bad_pixel_mask_value, + self.processed_data_gpu, + ), + ) + + def load_constant(self, constant_type, constant_data): + if constant_type is LpdConstants.Offset: + self.offset_map_gpu.set( + np.transpose( + constant_data.astype(np.float32), + (2, 0, 1, 3), + ) + ) + elif constant_type in (LpdConstants.BadPixelsDark, LpdConstants.BadPixelsFF): + self.bad_pixel_map_gpu |= cupy.asarray( + np.transpose( + constant_data.astype(np.uint32), + (2, 0, 1, 3), + ) + ) + elif constant_type is LpdConstants.GainAmpMap: + self.gain_amp_map_gpu.set( + np.transpose( + constant_data.astype(np.float32), + (2, 0, 1, 3), + ) + ) + elif constant_type is LpdConstants.FFMap: + self.gain_amp_map_gpu.set( + np.transpose( + constant_data.astype(np.float32), + (2, 0, 1, 3), + ) + ) + elif constant_type is LpdConstants.RelativeGain: + # TODO: figure out where to get + self.rel_gain_slopes_map_gpu.set( + np.transpose( + constant_data.astype(np.float32), + (2, 0, 1, 3), + ) + ) + """ + self.rel_gain_intercept_map_gpu.set( + np.transpose( + constant_data.astype(np.float32), + (2, 0, 1, 3), + ) + ) + """ + + def _init_kernels(self): + kernel_source = self._kernel_template.render( + { + "pixels_x": self.pixels_x, + "pixels_y": self.pixels_y, + "data_memory_cells": self.memory_cells, + "constant_memory_cells": self.constant_memory_cells, + "input_data_dtype": utils.np_dtype_to_c_type(self.input_data_dtype), + "output_data_dtype": utils.np_dtype_to_c_type(self.output_data_dtype), + "corr_enum": utils.enum_to_c_template(CorrectionFlags), + } + ) + self.source_module = cupy.RawModule(code=kernel_source) + self.correction_kernel = self.source_module.get_function("correct") + + def flush_buffers(self): + self.offset_map_gpu.fill(0) + self.gain_amp_map_gpu.fill(1) + self.rel_gain_intercept_map_gpu.fill(0) + self.rel_gain_slopes_map_gpu.fill(1) + self.flatfield_map_gpu.fill(1) + + +class LpdCalcatFriend(base_calcat.BaseCalcatFriend): + _constant_enum_class = LpdConstants + + def __init__(self, device, *args, **kwargs): + super().__init__(device, *args, **kwargs) + self._constants_need_conditions = { + LpdConstants.Offset: self.dark_condition, + LpdConstants.BadPixelsDark: self.dark_condition, + LpdConstants.GainAmpMap: self.category_condition, + LpdConstants.FFMap: self.category_condition, + LpdConstants.RelativeGain: self.category_condition, + LpdConstants.BadPixelsFF: self.category_condition, + } + + @staticmethod + def add_schema( + schema, + managed_keys, + param_prefix="constantParameters", + status_prefix="foundConstants", + ): + super(LpdCalcatFriend, LpdCalcatFriend).add_schema( + schema, managed_keys, "LPD-Type", param_prefix, status_prefix + ) + + ( + OVERWRITE_ELEMENT(schema) + .key(f"{param_prefix}.pixelsX") + .setNewDefaultValue(256) + .commit(), + + OVERWRITE_ELEMENT(schema) + .key(f"{param_prefix}.pixelsY") + .setNewDefaultValue(256) + .commit(), + + OVERWRITE_ELEMENT(schema) + .key(f"{param_prefix}.memoryCells") + .setNewDefaultValue(512) + .commit(), + + OVERWRITE_ELEMENT(schema) + .key(f"{param_prefix}.biasVoltage") + .setNewDefaultValue(250) + .commit(), + ) + + ( + DOUBLE_ELEMENT(schema) + .key(f"{param_prefix}.feedbackCapacitor") + .assignmentOptional() + .defaultValue(5) + .reconfigurable() + .commit(), + + DOUBLE_ELEMENT(schema) + .key(f"{param_prefix}.photonEnergy") + .assignmentOptional() + .defaultValue(12) + .reconfigurable() + .commit(), + + DOUBLE_ELEMENT(schema) + .key(f"{param_prefix}.category") + .displayedName("Some category") + .description( + "I don't know what this exactly is supposed to mean. Looking at " + "existing constants, I saw a GainAmpMap with 'category 1.0' and " + "'Source Energy 9.2 keV' and another one with 'category 0.0' and " + "'Source Energy 12.0 keV'. So maybe they're related." + ) + .assignmentOptional() + .defaultValue(1) + .reconfigurable() + .commit(), + ) + managed_keys.add(f"{param_prefix}.feedbackCapacitor") + managed_keys.add(f"{param_prefix}.photonEnergy") + managed_keys.add(f"{param_prefix}.category") + + base_calcat.add_status_schema_from_enum(schema, status_prefix, LpdConstants) + + def dark_condition(self): + res = base_calcat.OperatingConditions() + res["Memory cells"] = self._get_param("memoryCells") + res["Sensor Bias Voltage"] = self._get_param("biasVoltage") + res["Pixels X"] = self._get_param("pixelsX") + res["Pixels Y"] = self._get_param("pixelsY") + res["Feedback capacitor"] = self._get_param("feedbackCapacitor") + + return res + + def illuminated_condition(self): + res = self.dark_condition() + res["Source Energy"] = self._get_param("photonEnergy") + return res + + def category_condition(self): + res = self.illuminated_condition() + res["category"] = self._get_param("category") + return res + + +@KARABO_CLASSINFO("LpdCorrection", deviceVersion) +class LpdCorrection(BaseCorrection): + _correction_flag_class = CorrectionFlags + _correction_field_names = ( + ("offset", CorrectionFlags.OFFSET), + ("gainAmp", CorrectionFlags.GAIN_AMP), + ("relGain", CorrectionFlags.REL_GAIN), + ("flatfield", CorrectionFlags.FF_CORR), + ("badPixels", CorrectionFlags.BPMASK), + ) + _kernel_runner_class = LpdGpuRunner + _calcat_friend_class = LpdCalcatFriend + _constant_enum_class = LpdConstants + _managed_keys = BaseCorrection._managed_keys.copy() + + @staticmethod + def expectedParameters(expected): + ( + OVERWRITE_ELEMENT(expected) + .key("dataFormat.pixelsX") + .setNewDefaultValue(256) + .commit(), + + OVERWRITE_ELEMENT(expected) + .key("dataFormat.pixelsY") + .setNewDefaultValue(256) + .commit(), + + OVERWRITE_ELEMENT(expected) + .key("dataFormat.memoryCells") + .setNewDefaultValue(512) + .commit(), + + # TODO: determine ideal default + OVERWRITE_ELEMENT(expected) + .key("preview.selectionMode") + .setNewDefaultValue("frame") + .commit(), + ) + + LpdCalcatFriend.add_schema(expected, LpdCorrection._managed_keys) + add_correction_step_schema( + expected, LpdCorrection._managed_keys, LpdCorrection._correction_field_names + ) + + # mandatory: manager needs this in schema + ( + VECTOR_STRING_ELEMENT(expected) + .key("managedKeys") + .assignmentOptional() + .defaultValue(list(LpdCorrection._managed_keys)) + .commit() + ) + + @property + def input_data_shape(self): + return ( + self.unsafe_get("dataFormat.memoryCells"), + 1, + self.unsafe_get("dataFormat.pixelsY"), + self.unsafe_get("dataFormat.pixelsX"), + ) + + def _load_constant_to_runner(self, constant, constant_data): + self.kernel_runner.load_constant(constant, constant_data) + + def process_data( + self, + data_hash, + metadata, + source, + train_id, + image_data, + cell_table, + do_generate_preview, + ): + pulse_table = np.ravel(data_hash.get("image.pulseId")) + if self._frame_filter is not None: + try: + cell_table = cell_table[self._frame_filter] + pulse_table = pulse_table[self._frame_filter] + image_data = image_data[self._frame_filter] + except IndexError: + self.log_status_warn( + "Failed to apply frame filter, please check that it is valid!" + ) + return + + try: + # drop the singleton dimension for kernel runner + self.kernel_runner.load_data(image_data) + except ValueError as e: + self.log_status_warn(f"Failed to load data: {e}") + return + except Exception as e: + self.log_status_warn(f"Unknown exception when loading data to GPU: {e}") + + buffer_handle, buffer_array = self._shmem_buffer.next_slot() + self.kernel_runner.load_cell_table(cell_table) + self.kernel_runner.correct(self._correction_flag_enabled) + self.kernel_runner.reshape( + output_order=self.unsafe_get("dataFormat.outputAxisOrder"), + out=buffer_array, + ) + if do_generate_preview: + if self._correction_flag_enabled != self._correction_flag_preview: + self.kernel_runner.correct(self._correction_flag_preview) + ( + preview_slice_index, + preview_cell, + preview_pulse, + ) = utils.pick_frame_index( + self.unsafe_get("preview.selectionMode"), + self.unsafe_get("preview.index"), + cell_table, + pulse_table, + warn_func=self.log_status_warn, + ) + preview_raw, preview_corrected = self.kernel_runner.compute_previews( + preview_slice_index, + ) + + data_hash.set(self._image_data_path, buffer_handle) + data_hash.set(self._cell_table_path, cell_table[:, np.newaxis]) + data_hash.set("image.pulseId", pulse_table[:, np.newaxis]) + data_hash.set("calngShmemPaths", [self._image_data_path]) + self._write_output(data_hash, metadata) + if do_generate_preview: + self._write_preview_outputs( + ( + ("preview.outputRaw", preview_raw), + ("preview.outputCorrected", preview_corrected), + ), + metadata, + ) diff --git a/src/calng/kernels/lpd_gpu.cu b/src/calng/kernels/lpd_gpu.cu new file mode 100644 index 0000000000000000000000000000000000000000..3f2a20867e2eead43570993610479df707a8b342 --- /dev/null +++ b/src/calng/kernels/lpd_gpu.cu @@ -0,0 +1,75 @@ +#include <cuda_fp16.h> + +{{corr_enum}} + +extern "C" { + __global__ void correct(const {{input_data_dtype}}* data, // shape: memory cell, 1, y, x (I think) + const unsigned short* cell_table, + const unsigned char corr_flags, + const float* offset_map, // shape: cell, y, x, gain + const float* gain_amp_map, + const float* rel_gain_intercept_map, + const float* rel_gain_slopes_map, + const float* flatfield_map, + const unsigned int* bad_pixel_map, + const float bad_pixel_mask_value, + {{output_data_dtype}}* output) { + const size_t X = {{pixels_x}}; + const size_t Y = {{pixels_y}}; + const size_t memory_cells = {{data_memory_cells}}; + const size_t map_memory_cells = {{constant_memory_cells}}; + + const size_t memory_cell = blockIdx.x * blockDim.x + threadIdx.x; + const size_t y = blockIdx.y * blockDim.y + threadIdx.y; + const size_t x = blockIdx.z * blockDim.z + threadIdx.z; + + if (memory_cell >= memory_cells || y >= Y || x >= X) { + return; + } + + const size_t data_stride_x = 1; + const size_t data_stride_y = X * data_stride_x; + const size_t data_stride_cell = Y * data_stride_y; + const size_t data_index = memory_cell * data_stride_cell + y * data_stride_y + x * data_stride_x; + const unsigned short raw_data_value = data[data_index]; + const unsigned char gain = raw_data_value >> 12; + float corrected = (float)(raw_data_value & 0xfff); + + const size_t gm_map_stride_gain = 1; + const size_t gm_map_stride_x = 3 * gm_map_stride_gain; + const size_t gm_map_stride_y = X * gm_map_stride_x; + const size_t gm_map_stride_cell = Y * gm_map_stride_y; + + const size_t map_cell = cell_table[memory_cell]; + if (map_cell < map_memory_cells) { + const size_t gm_map_index = gain * gm_map_stride_gain + + map_cell * gm_map_stride_cell + + y * gm_map_stride_y + + x * gm_map_stride_x; + + if ((corr_flags & BPMASK) && bad_pixel_map[gm_map_index]) { + corrected = bad_pixel_mask_value; + } else { + // TODO: save GPU memory by combining maps + if (corr_flags & OFFSET) { + corrected -= offset_map[gm_map_index]; + } + if (corr_flags & GAIN_AMP) { + corrected *= gain_amp_map[gm_map_index]; + } + if (corr_flags & REL_GAIN) { + corrected -= rel_gain_intercept_map[gm_map_index]; + corrected *= rel_gain_slopes_map[gm_map_index]; + } + if (corr_flags & FF_CORR) { + corrected *= flatfield_map[gm_map_index]; + } + } + } + {% if output_data_dtype == "half" %} + output[data_index] = __float2half(corrected); + {% else %} + output[data_index] = ({{output_data_dtype}})corrected; + {% endif %} + } +} diff --git a/src/calng/mdl_geometry_base.py b/src/calng/mdl_geometry_base.py index c56125f733a6caee568811d437c6cd2fa9fbaf32..4115e6038263295886a40f498d86c2f3b1724b0c 100644 --- a/src/calng/mdl_geometry_base.py +++ b/src/calng/mdl_geometry_base.py @@ -260,11 +260,35 @@ class MdlAgipd1MGeometry(ManualQuadrantsGeometryBase): ) +class MdlDssc1MGeometry(ManualQuadrantsGeometryBase): + geometry_class = extra_geom.DSSC_1MGeometry + quadrantCorners = makeQuadrantCornersNode( + [ + (-130, 5), + (-130, -125), + (5, -125), + (5, 5), + ] + ) + + +class MdlLpd1MGeometry(ManualQuadrantsGeometryBase): + geometry_class = extra_geom.LPD_1MGeometry + quadrantCorners = makeQuadrantCornersNode( + [ + (11.4, 299), + (-11.5, 8), + (254.5, -16), + (278.5, 275), + ] + ) + + class ModuleListItem(Configurable): - posX = Double(defaultValue=0) - posY = Double(defaultValue=0) - orientX = Int32(defaultValue=1) - orientY = Int32(defaultValue=1) + posX = Double() + posY = Double() + orientX = Int32() + orientY = Int32() class ManualModuleListGeometryBase(ManualGeometryBase): diff --git a/src/calng/utils.py b/src/calng/utils.py index 0e560b053536736a6663c4bb44d313229371f4f0..b14e36b84675cd06a62824828b158edcd9edbdd4 100644 --- a/src/calng/utils.py +++ b/src/calng/utils.py @@ -1,3 +1,4 @@ +import enum import collections import functools import inspect @@ -347,3 +348,31 @@ class ChainHash: if h.has(key): return h[key] raise KeyError() + + +class BadPixelValues(enum.IntFlag): + """The European XFEL Bad Pixel Encoding + + Straight from pycalibration's enum.py""" + + OFFSET_OUT_OF_THRESHOLD = 2 ** 0 + NOISE_OUT_OF_THRESHOLD = 2 ** 1 + OFFSET_NOISE_EVAL_ERROR = 2 ** 2 + NO_DARK_DATA = 2 ** 3 + CI_GAIN_OUT_OF_THRESHOLD = 2 ** 4 + CI_LINEAR_DEVIATION = 2 ** 5 + CI_EVAL_ERROR = 2 ** 6 + FF_GAIN_EVAL_ERROR = 2 ** 7 + FF_GAIN_DEVIATION = 2 ** 8 + FF_NO_ENTRIES = 2 ** 9 + CI2_EVAL_ERROR = 2 ** 10 + VALUE_IS_NAN = 2 ** 11 + VALUE_OUT_OF_RANGE = 2 ** 12 + GAIN_THRESHOLDING_ERROR = 2 ** 13 + DATA_STD_IS_ZERO = 2 ** 14 + ASIC_STD_BELOW_NOISE = 2 ** 15 + INTERPOLATED = 2 ** 16 + NOISY_ADC = 2 ** 17 + OVERSCAN = 2 ** 18 + NON_SENSITIVE = 2 ** 19 + NON_LIN_RESPONSE_REGION = 2 ** 20