From 19934f5ac68c371b51e60dbd1ec094643093c44c Mon Sep 17 00:00:00 2001 From: David Hammer <dhammer@mailbox.org> Date: Fri, 30 Jun 2023 11:25:59 +0200 Subject: [PATCH] Squashed commit of the following: commit aaa8a45d79a7aeb0a34bc48db1b1f1b423953133 Merge: 8033df8 c838797 Author: David Hammer <david.hammer@xfel.eu> Date: Fri Jun 30 09:24:26 2023 +0000 Merge branch 'master' into 'pnccd-support' # Conflicts: # setup.py # src/calng/base_geometry.py # src/calng/geometries/__init__.py commit 8033df8a239a3c7ec63b34831b101f9339bdccb5 Author: David Hammer <dhammer@mailbox.org> Date: Fri Jun 30 11:19:06 2023 +0200 Fix scene layout commit cb2951a8b90b8bce45fac39064188c9ac20778ef Author: David Hammer <dhammer@mailbox.org> Date: Fri Jun 30 11:13:32 2023 +0200 Prefer pass over ellipsis for empty body commit 4eb78216789fb125f983ca57ddaa11a044b44da4 Author: David Hammer <dhammer@mailbox.org> Date: Tue Jun 27 15:59:21 2023 +0200 Use faster common mode (row / col) for ePix100, too commit df7593edad4be702be5cddab3b735d787a951844 Author: David Hammer <dhammer@mailbox.org> Date: Tue Jun 27 15:39:40 2023 +0200 Tweak custom common mode implementation, add test commit fcc6a0899004bf6b8a985ceb8f4a323e5b644513 Author: David Hammer <dhammer@mailbox.org> Date: Tue Jun 27 14:07:14 2023 +0200 Fix correct/reshape/correct for ePix100, too commit a35ea8173b4e7c3463cb6c3ccde7ebed2b0f63db Author: David Hammer <dhammer@mailbox.org> Date: Tue Jun 27 14:04:13 2023 +0200 Fix order of correct, reshape, correct commit bb89405a14bd74e41af36bbf81f9d583b741531b Author: David Hammer <dhammer@mailbox.org> Date: Tue Jun 27 10:28:41 2023 +0200 Don't always claim to load to GPU commit 156a07931af7a83c35ec804dd11118476b8f0b4a Author: David Hammer <dhammer@mailbox.org> Date: Mon Jun 26 18:13:03 2023 +0200 Fix miscellaneous issues, properly add geometry commit b773af39b533e69710e67f86486675df27fe4931 Author: David Hammer <dhammer@mailbox.org> Date: Mon Jun 26 13:47:05 2023 +0200 Add pnCCD geometry class, use Configurable.get_root instead of ad-hoc implementation commit e6c0399391f96df49e013e89069e0a78b9743c43 Author: David Hammer <dhammer@mailbox.org> Date: Mon Jun 26 13:46:39 2023 +0200 Fix typo in UINT64_ELEMENT commit 1b408d3011d598a35e62d4a2c0a7e5a61c27557d Author: David Hammer <dhammer@mailbox.org> Date: Thu Jun 1 14:43:49 2023 +0200 By default, don't ignore bias voltage commit bdecca450bc6756f19174d91e9574bc4960f46e3 Author: David Hammer <dhammer@mailbox.org> Date: Fri May 26 17:52:41 2023 +0200 Fix output schema commit b405d5a1356eb22e10c99d321570dcd3e6a7f65b Author: David Hammer <dhammer@mailbox.org> Date: Fri May 26 17:28:32 2023 +0200 Call it ignoreBiasVoltage instead commit 5c45ef92b3f8ea6ae4fd95f28e27f77b665c9078 Author: David Hammer <dhammer@mailbox.org> Date: Fri May 26 15:04:07 2023 +0200 Add faster Cython common mode commit c7c62e2b3b8475b1d93d2dbdcb9d0d9e4945b05e Author: sqsonc <sqsonc@exflong22.desy.de> Date: Tue May 23 10:43:48 2023 +0200 Put output data in proper hash keys commit 0779b7e537e6bf8d14c4159f7469584ec0ab79b7 Author: David Hammer <dhammer@mailbox.org> Date: Fri May 12 14:09:46 2023 +0200 Removing block mode CM for now commit f322ae74039b021757135727678c341dcb63f995 Author: David Hammer <dhammer@mailbox.org> Date: Fri May 12 12:51:45 2023 +0200 Add WIP pnCCD correction device --- setup.py | 11 +- src/calng/RoiTool.py | 3 +- src/calng/base_correction.py | 2 +- src/calng/base_geometry.py | 117 ++++- src/calng/corrections/Epix100Correction.py | 29 +- src/calng/corrections/PnccdCorrection.py | 493 +++++++++++++++++++++ src/calng/geometries/PnccdGeometry.py | 19 + src/calng/geometries/__init__.py | 2 + src/calng/kernels/common_cpu.pyx | 70 +++ src/calng/scenes.py | 82 +++- src/calng/schemas.py | 31 ++ src/tests/test_pnccd_kernels.py | 50 +++ 12 files changed, 881 insertions(+), 28 deletions(-) create mode 100644 src/calng/corrections/PnccdCorrection.py create mode 100644 src/calng/geometries/PnccdGeometry.py create mode 100644 src/calng/kernels/common_cpu.pyx create mode 100644 src/tests/test_pnccd_kernels.py diff --git a/setup.py b/setup.py index 2f94de91..92c44ccc 100644 --- a/setup.py +++ b/setup.py @@ -30,6 +30,7 @@ setup(name='calng', 'Epix100Correction = calng.corrections.Epix100Correction:Epix100Correction', 'Gotthard2Correction = calng.corrections.Gotthard2Correction:Gotthard2Correction', 'JungfrauCorrection = calng.corrections.JungfrauCorrection:JungfrauCorrection', + 'PnccdCorrection = calng.corrections.PnccdCorrection:PnccdCorrection', 'LpdCorrection = calng.corrections.LpdCorrection:LpdCorrection', 'LpdminiCorrection = calng.corrections.LpdminiCorrection:LpdminiCorrection', 'ShmemToZMQ = calng.ShmemToZMQ:ShmemToZMQ', @@ -48,9 +49,10 @@ setup(name='calng', 'Agipd500KGeometry = calng.geometries.Agipd500KGeometry:Agipd500KGeometry', 'Dssc1MGeometry = calng.geometries:Dssc1MGeometry.Dssc1MGeometry', 'Epix100Geometry = calng.geometries:Epix100Geometry.Epix100Geometry', + 'JungfrauGeometry = calng.geometries:JungfrauGeometry.JungfrauGeometry', 'Lpd1MGeometry = calng.geometries:Lpd1MGeometry.Lpd1MGeometry', 'LpdminiGeometry = calng.geometries:LpdminiGeometry.LpdminiGeometry', - 'JungfrauGeometry = calng.geometries:JungfrauGeometry.JungfrauGeometry', + 'PnccdGeometry = calng.geometries:PnccdGeometry.PnccdGeometry', 'RoiTool = calng.RoiTool:RoiTool', ], }, @@ -81,6 +83,13 @@ setup(name='calng', extra_compile_args=['-O3', '-march=native', '-fopenmp'], extra_link_args=['-fopenmp'], ), + Extension( + 'calng.kernels.common_cython', + ['src/calng/kernels/common_cpu.pyx'], + extra_compile_args=['-O3', '-march=native', '-fopenmp'], + extra_link_args=['-fopenmp'], + language="c++", + ), Extension( 'calng.kernels.lpd_cython', ['src/calng/kernels/lpd_cpu.pyx'], diff --git a/src/calng/RoiTool.py b/src/calng/RoiTool.py index 4618aebe..4f772c96 100644 --- a/src/calng/RoiTool.py +++ b/src/calng/RoiTool.py @@ -54,8 +54,7 @@ class PreviewSettingsNode(Configurable): @Double(accessMode=AccessMode.RECONFIGURABLE, defaultValue=2, unitSymbol=Unit.HERTZ) def maxPreviewRate(self, value): self.maxPreviewRate = value - parent = self.get_root() - parent._throttler = utils.SkippingThrottler(1 / self.maxPreviewRate.value) + self.get_root()._throttler = utils.SkippingThrottler(1 / self.maxPreviewRate.value) class HistogramSettingsNode(Configurable): diff --git a/src/calng/base_correction.py b/src/calng/base_correction.py index 57a62cee..a205471d 100644 --- a/src/calng/base_correction.py +++ b/src/calng/base_correction.py @@ -102,7 +102,7 @@ class BaseCorrection(PythonDevice): self.set(key, True) self._update_correction_flags() - self.log_status_info(f"Done loading {constant.name} to GPU") + self.log_status_info(f"Done loading {constant.name} to runner") @property def input_data_shape(self): diff --git a/src/calng/base_geometry.py b/src/calng/base_geometry.py index 5ed201da..f7c70e38 100644 --- a/src/calng/base_geometry.py +++ b/src/calng/base_geometry.py @@ -1,4 +1,5 @@ import contextlib +import functools import logging import matplotlib.pyplot as plt @@ -54,18 +55,66 @@ def make_x_y_coordinate_node( return Node(XYCoordinate, **({} if node_args is None else node_args)) -def make_x_y_offset_node(): - return make_x_y_coordinate_node( +def make_x_y_z_coordinate_node( + default_x, + default_y, + default_z, + x_args=None, + y_args=None, + z_args=None, + node_args=None, +): + class XYZCoordinate(Configurable): + x = Double( + defaultValue=default_x, + accessMode=AccessMode.RECONFIGURABLE, + assignment=Assignment.OPTIONAL, + **({} if x_args is None else x_args), + ) + y = Double( + defaultValue=default_y, + accessMode=AccessMode.RECONFIGURABLE, + assignment=Assignment.OPTIONAL, + **({} if y_args is None else y_args), + ) + z = Double( + defaultValue=default_z, + accessMode=AccessMode.RECONFIGURABLE, + assignment=Assignment.OPTIONAL, + **({} if z_args is None else z_args), + ) + + return Node(XYZCoordinate, **({} if node_args is None else node_args)) + + +make_x_y_offset_node = functools.partial( + make_x_y_coordinate_node, + 0, + 0, + x_args={"unitSymbol": Unit.METER}, + y_args={"unitSymbol": Unit.METER}, + node_args={ + "displayedName": "Offset", + "description": "See EXtra-geom documentation for details. This offset is " + "applied to entire detector after initial geometry is created from manual " + "parameters. Example: To move entire geometry up by 2 mm relative to beam, " + "set offset.y to 2e-3.", + }, +) + + +def make_x_y_z_relative_offset_node(title): + return make_x_y_z_coordinate_node( + 0, 0, 0, x_args={"unitSymbol": Unit.METER}, y_args={"unitSymbol": Unit.METER}, + z_args={"unitSymbol": Unit.METER}, node_args={ - "displayedName": "Offset", - "description": "See EXtra-geom documentation for details. This offset is " - "applied to entire detector after initial geometry is created from manual " - "parameters. Example: To move entire geometry up by 2 mm relative to beam, " - "set offset.y to 2e-3.", + "displayedName": title, + "description": "See EXtra-geom documentation for details. This is used to " + "set the relative offset for the two halves of a pnCCD", }, ) @@ -383,6 +432,60 @@ class ManualOriginGeometryBase(ManualGeometryBase): raise NotImplementedError() +class RelativePositionNode(BaseManualGeometryConfigNode): + gap = Double( + defaultValue=0.004, + accessMode=AccessMode.RECONFIGURABLE, + assignment=Assignment.OPTIONAL, + unitSymbol=Unit.METER, + displayedName="Gap", + description="Gap between modules. 4mm (~50) by default.", + ) + topOffset = make_x_y_z_relative_offset_node("Top offset") + bottomOffset = make_x_y_z_relative_offset_node("Bottom offset") + + +class ManualRelativePositionsGeometryBase(ManualGeometryBase): + """Only used for pnCCD for now: specify geometry primarily via the gap between the + two detector modules. Optionally tweak by offsetting top and / or bottom half. If + this is useful for other detectors at some point, should refactor to set individual + defaults.""" + + manualSetting = Node(RelativePositionNode) + + async def _set_from_manual_config(self): + self._set_status("Updating geometry from manual configuration") + geometry = self.geometry_class.from_relative_positions( + gap=self.manualSetting.gap.value, + top_offset=( + self.manualSetting.topOffset.x.value, + self.manualSetting.topOffset.y.value, + self.manualSetting.topOffset.z.value, + ), + bottom_offset=( + self.manualSetting.bottomOffset.x.value, + self.manualSetting.bottomOffset.y.value, + self.manualSetting.bottomOffset.z.value, + ), + ) + await self._set_geometry(geometry) + + @slot + def requestScene(self, params): + name = params.get("name", default="overview") + if name == "overview": + scene_data = scenes.relative_geometry_overview( + self.deviceId, + self.getDeviceSchema(), + ) + payload = Hash("success", True, "name", name, "data", scene_data) + + return Hash("type", "deviceScene", "origin", self.deviceId, "payload", payload) + + async def _update_manual_from_current(self): + raise NotImplementedError() + + def make_quadrant_corners_node(default_values): assert len(default_values) == 4 assert all(len(x) == 2 for x in default_values) diff --git a/src/calng/corrections/Epix100Correction.py b/src/calng/corrections/Epix100Correction.py index 3d72b18d..6e16198b 100644 --- a/src/calng/corrections/Epix100Correction.py +++ b/src/calng/corrections/Epix100Correction.py @@ -21,6 +21,7 @@ from .. import ( utils, ) from .._version import version as deviceVersion +from ..kernels import common_cython class Constants(enum.Enum): @@ -251,11 +252,11 @@ class Epix100CpuRunner(base_kernel_runner.BaseKernelRunner): if flags & CorrectionFlags.COMMONMODE: # per rectangular block that looks like something is going on - masked = np.ma.masked_array( - data=output, - mask=(self._q_bad_pixel_map[q] != 0) - | (output > self._q_noise_map[q] * cm_noise_sigma), - ) + cm_mask=( + (self._q_bad_pixel_map[q] != 0) + | (output > self._q_noise_map[q] * cm_noise_sigma) + ).astype(np.uint8, copy=False) + masked = np.ma.masked_array(data=output, mask=cm_mask) if cm_block: for block in np.hsplit(masked, 4): if block.count() < block.size * cm_min_frac: @@ -263,14 +264,10 @@ class Epix100CpuRunner(base_kernel_runner.BaseKernelRunner): block.data[:] -= np.ma.median(block) if cm_row: - subset_rows = masked.count(axis=1) >= masked.shape[1] * cm_min_frac - output[subset_rows] -= np.ma.median( - masked[subset_rows], axis=1, keepdims=True - ) + common_cython.cm_fs(output, cm_mask, cm_noise_sigma, cm_min_frac) if cm_col: - subset_cols = masked.count(axis=0) >= masked.shape[0] * cm_min_frac - output[:, subset_cols] -= np.ma.median(masked[:, subset_cols], axis=0) + common_cython.cm_ss(output, cm_mask, cm_noise_sigma, cm_min_frac) if flags & CorrectionFlags.RELGAIN: output *= self._q_rel_gain_map[q] @@ -302,7 +299,7 @@ class Epix100CpuRunner(base_kernel_runner.BaseKernelRunner): ), range(4), ): - ... + pass # just run through to await map @KARABO_CLASSINFO("Epix100Correction", deviceVersion) @@ -454,15 +451,15 @@ class Epix100Correction(base_correction.BaseCorrection): self.kernel_runner.correct( flags=self._correction_flag_enabled, **args_which_should_be_cached ) + self.kernel_runner.reshape( + output_order=self.unsafe_get("dataFormat.outputAxisOrder"), + out=buffer_array, + ) if self._correction_flag_enabled != self._correction_flag_preview: self.kernel_runner.correct( flags=self._correction_flag_preview, **args_which_should_be_cached, ) - self.kernel_runner.reshape( - output_order=self.unsafe_get("dataFormat.outputAxisOrder"), - out=buffer_array, - ) if self._use_shmem_handles: # TODO: consider better name for key... diff --git a/src/calng/corrections/PnccdCorrection.py b/src/calng/corrections/PnccdCorrection.py new file mode 100644 index 00000000..b2adb55d --- /dev/null +++ b/src/calng/corrections/PnccdCorrection.py @@ -0,0 +1,493 @@ +import concurrent.futures +import enum +import functools + +import numpy as np +from karabo.bound import ( + BOOL_ELEMENT, + DOUBLE_ELEMENT, + KARABO_CLASSINFO, + OUTPUT_CHANNEL, + OVERWRITE_ELEMENT, + VECTOR_STRING_ELEMENT, + Schema, +) + +from .. import ( + base_calcat, + base_correction, + base_kernel_runner, + schemas, + utils, +) +from .._version import version as deviceVersion +from ..kernels import common_cython + + +class Constants(enum.Enum): + BadPixelsDarkCCD = enum.auto() + OffsetCCD = enum.auto() + NoiseCCD = enum.auto() + RelativeGainCCD = enum.auto() + CTICCD = enum.auto() # what is this? + + +class CorrectionFlags(enum.IntFlag): + NONE = 0 + OFFSET = 2**1 + COMMONMODE = 2**2 + RELGAIN = 2**3 + BPMASK = 2**4 + + +class PnccdCalcatFriend(base_calcat.BaseCalcatFriend): + _constant_enum_class = Constants + + @property + def _constants_need_conditions(self): + return { + Constants.OffsetCCD: self.dark_condition, + Constants.BadPixelsDarkCCD: self.dark_condition, + Constants.NoiseCCD: self.dark_condition, + Constants.RelativeGainCCD: self.illuminated_condition, + } + + @staticmethod + def add_schema(schema, managed_keys): + super(PnccdCalcatFriend, PnccdCalcatFriend).add_schema( + schema, managed_keys, "pnCCD-Type" + ) + + # set some defaults for common parameters + ( + OVERWRITE_ELEMENT(schema) + .key("constantParameters.memoryCells") + .setNewDefaultValue(1) + .commit(), + + OVERWRITE_ELEMENT(schema) + .key("constantParameters.biasVoltage") + .setNewDefaultValue(300) + .commit(), + + OVERWRITE_ELEMENT(schema) + .key("constantParameters.pixelsX") + .setNewDefaultValue(1024) + .commit(), + + OVERWRITE_ELEMENT(schema) + .key("constantParameters.pixelsY") + .setNewDefaultValue(1024) + .commit(), + ) + + ( + BOOL_ELEMENT(schema) + .key("constantParameters.ignoreBiasVoltage") + .displayedName("Ignore bias voltage parameter") + .assignmentOptional() + .defaultValue(False) + .reconfigurable() + .commit(), + + DOUBLE_ELEMENT(schema) + .key("constantParameters.integrationTime") + .displayedName("Integration Time") + .assignmentOptional() + .defaultValue(70) + .reconfigurable() + .commit(), + + DOUBLE_ELEMENT(schema) + .key("constantParameters.sensorTemperature") + .assignmentOptional() + .defaultValue(233) + .reconfigurable() + .commit(), + + DOUBLE_ELEMENT(schema) + .key("constantParameters.gainSetting") + .displayedName("Gain Setting") + .assignmentOptional() + .defaultValue(16) + .reconfigurable() + .commit(), + + DOUBLE_ELEMENT(schema) + .key("constantParameters.sourceEnergy") + .displayedName("Source Energy") + .assignmentOptional() + .defaultValue(1.6) + .reconfigurable() + .commit(), + ) + + managed_keys.add("constantParameters.ignoreBiasVoltage") + managed_keys.add("constantParameters.integrationTime") + managed_keys.add("constantParameters.gainSetting") + managed_keys.add("constantParameters.sourceEnergy") + managed_keys.add("constantParameters.sensorTemperature") + + base_calcat.add_status_schema_from_enum(schema, Constants) + + def dark_condition(self): + res = base_calcat.OperatingConditions() + res["Memory cells"] = self._get_param("memoryCells") + if not self._get_param("ignoreBiasVoltage"): + res["Sensor Bias Voltage"] = self._get_param("biasVoltage") + res["Pixels X"] = self._get_param("pixelsX") + res["Pixels Y"] = self._get_param("pixelsY") + res["Integration Time"] = self._get_param("integrationTime") + res["Sensor Temperature"] = self._get_param("sensorTemperature") + res["Gain Setting"] = self._get_param("gainSetting") + + return res + + def illuminated_condition(self): + res = self.dark_condition() + + res["Source Energy"] = self._get_param("sourceEnergy") + + return res + + +class PnccdCpuRunner(base_kernel_runner.BaseKernelRunner): + _corrected_axis_order = "fxy" + _xp = np + _gpu_based = False + + @property + def input_shape(self): + return (self.pixels_x, self.pixels_y) + + @property + def preview_shape(self): + return (self.pixels_x, self.pixels_y) + + @property + def processed_shape(self): + return self.input_shape + + @property + def map_shape(self): + return (self.pixels_x, self.pixels_y) + + def __init__( + self, + pixels_x, + pixels_y, + frames, # will be 1, will be ignored + constant_memory_cells, + input_data_dtype=np.uint16, + output_data_dtype=np.float32, + ): + assert ( + output_data_dtype == np.float32 + ), "Alternative output types not supported yet" + + super().__init__( + pixels_x, + pixels_y, + frames, + constant_memory_cells, + output_data_dtype, + ) + + self.input_data = None + self.processed_data = np.empty(self.processed_shape, dtype=np.float32) + + self.offset_map = np.empty(self.map_shape, dtype=np.float32) + self.rel_gain_map = np.empty(self.map_shape, dtype=np.float32) + self.bad_pixel_map = np.empty(self.map_shape, dtype=np.uint32) + self.noise_map = np.empty(self.map_shape, dtype=np.float32) + # will do everything by quadrant + self.thread_pool = concurrent.futures.ThreadPoolExecutor(max_workers=4) + self._q_input_data = None + self._q_processed_data = utils.quadrant_views(self.processed_data) + self._q_offset_map = utils.quadrant_views(self.offset_map) + self._q_rel_gain_map = utils.quadrant_views(self.rel_gain_map) + self._q_bad_pixel_map = utils.quadrant_views(self.bad_pixel_map) + self._q_noise_map = utils.quadrant_views(self.noise_map) + + def __del__(self): + self.thread_pool.shutdown() + + @property + def preview_data_views(self): + return (self.input_data, self.processed_data) + + def load_data(self, image_data): + # should almost never squeeze, but ePix doesn't do burst mode, right? + self.input_data = image_data.astype(np.uint16, copy=False).squeeze() + self._q_input_data = utils.quadrant_views(self.input_data) + + def load_constant(self, constant_type, data): + if constant_type is Constants.OffsetCCD: + self.offset_map[:] = data.squeeze().astype(np.float32) + elif constant_type is Constants.RelativeGainCCD: + self.rel_gain_map[:] = data.squeeze().astype(np.float32) + elif constant_type is Constants.BadPixelsDarkCCD: + self.bad_pixel_map[:] = data.squeeze() + elif constant_type is Constants.NoiseCCD: + self.noise_map[:] = data.squeeze() + else: + raise Exception(f"Unknown constant type {constant_type}") + + def flush_buffers(self, constants): + if Constants.OffsetCCD in constants: + self.offset_map.fill(0) + if Constants.RelativeGainCCD in constants: + self.rel_gain_map.fill(1) + if constants & {Constants.BadPixelsDarkCCD}: + self.bad_pixel_map.fill(0) + if Constants.NoiseCCD in constants: + self.noise_map.fill(np.inf) + + def _correct_quadrant( + self, + q, + flags, + bad_pixel_mask_value, + cm_noise_sigma, + cm_min_frac, + cm_row, + cm_col, + ): + output = self._q_processed_data[q] + output[:] = self._q_input_data[q].astype(np.float32) + if flags & CorrectionFlags.OFFSET: + output -= self._q_offset_map[q] + + if flags & CorrectionFlags.COMMONMODE: + cm_mask = ( + (self._q_bad_pixel_map[q] != 0) + | (output > self._q_noise_map[q] * cm_noise_sigma) + ).astype(np.uint8, copy=False) + if cm_row: + common_cython.cm_fs( + output, + cm_mask, + cm_noise_sigma, + cm_min_frac, + ) + + if cm_col: + common_cython.cm_ss( + output, + cm_mask, + cm_noise_sigma, + cm_min_frac, + ) + + if flags & CorrectionFlags.RELGAIN: + output *= self._q_rel_gain_map[q] + + if flags & CorrectionFlags.BPMASK: + output[self._q_bad_pixel_map[q] != 0] = bad_pixel_mask_value + + def correct( + self, + flags, + bad_pixel_mask_value=np.nan, + cm_noise_sigma=5, + cm_min_frac=0.25, + cm_row=True, + cm_col=True, + ): + # NOTE: how to best clean up all these duplicated parameters? + for result in self.thread_pool.map( + functools.partial( + self._correct_quadrant, + flags=flags, + bad_pixel_mask_value=bad_pixel_mask_value, + cm_noise_sigma=cm_noise_sigma, + cm_min_frac=cm_min_frac, + cm_row=cm_row, + cm_col=cm_col, + ), + range(4), + ): + pass # just run through to await map + + +@KARABO_CLASSINFO("PnccdCorrection", deviceVersion) +class PnccdCorrection(base_correction.BaseCorrection): + _correction_flag_class = CorrectionFlags + _correction_steps = ( + ("offset", CorrectionFlags.OFFSET, {Constants.OffsetCCD}), + ("relGain", CorrectionFlags.RELGAIN, {Constants.RelativeGainCCD}), + ("commonMode", CorrectionFlags.COMMONMODE, {Constants.NoiseCCD}), + ("badPixels", CorrectionFlags.BPMASK, {Constants.BadPixelsDarkCCD}), + ) + _image_data_path = "data.image" + _kernel_runner_class = PnccdCpuRunner + _calcat_friend_class = PnccdCalcatFriend + _constant_enum_class = Constants + _managed_keys = base_correction.BaseCorrection._managed_keys.copy() + _preview_outputs = ["outputRaw", "outputCorrected"] + _cell_table_path = None + _pulse_table_path = None + _warn_memory_cell_range = False + + @staticmethod + def expectedParameters(expected): + ( + OUTPUT_CHANNEL(expected) + .key("dataOutput") + .dataSchema(schemas.pnccd_output_schema(use_shmem_handle=False)) + .commit(), + + OVERWRITE_ELEMENT(expected) + .key("dataFormat.pixelsX") + .setNewDefaultValue(1024) + .commit(), + + OVERWRITE_ELEMENT(expected) + .key("dataFormat.pixelsY") + .setNewDefaultValue(1024) + .commit(), + + OVERWRITE_ELEMENT(expected) + .key("dataFormat.frames") + .setNewDefaultValue(1) + .commit(), + + # TODO: disable preview selection mode + + OVERWRITE_ELEMENT(expected) + .key("useShmemHandles") + .setNewDefaultValue(False) + .commit(), + + OVERWRITE_ELEMENT(expected) + .key("outputShmemBufferSize") + .setNewDefaultValue(2) + .commit(), + ) + + base_correction.add_preview_outputs( + expected, PnccdCorrection._preview_outputs + ) + base_correction.add_correction_step_schema( + expected, + PnccdCorrection._managed_keys, + PnccdCorrection._correction_steps, + ) + ( + DOUBLE_ELEMENT(expected) + .key("corrections.commonMode.noiseSigma") + .assignmentOptional() + .defaultValue(5) + .reconfigurable() + .commit(), + + DOUBLE_ELEMENT(expected) + .key("corrections.commonMode.minFrac") + .assignmentOptional() + .defaultValue(0.25) + .reconfigurable() + .commit(), + + BOOL_ELEMENT(expected) + .key("corrections.commonMode.enableRow") + .assignmentOptional() + .defaultValue(True) + .reconfigurable() + .commit(), + + BOOL_ELEMENT(expected) + .key("corrections.commonMode.enableCol") + .assignmentOptional() + .defaultValue(True) + .reconfigurable() + .commit(), + ) + PnccdCorrection._managed_keys |= { + "corrections.commonMode.noiseSigma", + "corrections.commonMode.minFrac", + "corrections.commonMode.enableRow", + "corrections.commonMode.enableCol", + } + PnccdCalcatFriend.add_schema(expected, PnccdCorrection._managed_keys) + # TODO: bad pixel node? + + # mandatory: manager needs this in schema + ( + VECTOR_STRING_ELEMENT(expected) + .key("managedKeys") + .assignmentOptional() + .defaultValue(list(PnccdCorrection._managed_keys)) + .commit() + ) + + @property + def input_data_shape(self): + # TODO: check + return ( + self.unsafe_get("dataFormat.pixelsX"), + self.unsafe_get("dataFormat.pixelsY"), + ) + + def __init__(self, config): + super().__init__(config) + if config.get("useShmemHandles", default=False): + def aux(): + schema_override = Schema() + ( + OUTPUT_CHANNEL(schema_override) + .key("dataOutput") + .dataSchema(schemas.pnccd_output_schema(use_shmem_handle=False)) + .commit(), + ) + self.updateSchema(schema_override) + + self.registerInitialFunction(aux) + + def process_data( + self, + data_hash, + metadata, + source, + train_id, + image_data, + cell_table, # will be None + pulse_table, # ditto + ): + self.kernel_runner.load_data(image_data) + + buffer_handle, buffer_array = self._shmem_buffer.next_slot() + args_which_should_be_cached = dict( + cm_noise_sigma=self.unsafe_get("corrections.commonMode.noiseSigma"), + cm_min_frac=self.unsafe_get("corrections.commonMode.minFrac"), + cm_row=self.unsafe_get("corrections.commonMode.enableRow"), + cm_col=self.unsafe_get("corrections.commonMode.enableCol"), + ) + self.kernel_runner.correct( + flags=self._correction_flag_enabled, **args_which_should_be_cached + ) + self.kernel_runner.reshape( + output_order=self.unsafe_get("dataFormat.outputAxisOrder"), + out=buffer_array, + ) + if self._correction_flag_enabled != self._correction_flag_preview: + self.kernel_runner.correct( + flags=self._correction_flag_preview, + **args_which_should_be_cached, + ) + + if self._use_shmem_handles: + data_hash.set(self._image_data_path, buffer_handle) + data_hash.set("calngShmemPaths", [self._image_data_path]) + else: + data_hash.set(self._image_data_path, buffer_array) + data_hash.set("calngShmemPaths", []) + + self._write_output(data_hash, metadata) + + # note: base class preview machinery assumes burst mode, shortcut it + self._preview_friend.write_outputs( + metadata, *self.kernel_runner.preview_data_views + ) + + def _load_constant_to_runner(self, constant, constant_data): + self.kernel_runner.load_constant(constant, constant_data) diff --git a/src/calng/geometries/PnccdGeometry.py b/src/calng/geometries/PnccdGeometry.py new file mode 100644 index 00000000..c3b84cca --- /dev/null +++ b/src/calng/geometries/PnccdGeometry.py @@ -0,0 +1,19 @@ +import extra_geom + +from ..base_geometry import ManualRelativePositionsGeometryBase + + +class FixedUpPnccdGeometry(extra_geom.PNCCDGeometry): + # DetectorAssembler neeeds expected_data_shape to match actual data shape + expected_data_shape = (1, 1024, 1024) + + # and these two methods need tweaking to avoid that breaking assertions + def _ensure_shape(self, data): + return data.reshape(self.expected_data_shape) + + def split_tiles(self, module_data): + return [module_data[:512], module_data[512:]] + + +class PnccdGeometry(ManualRelativePositionsGeometryBase): + geometry_class = FixedUpPnccdGeometry diff --git a/src/calng/geometries/__init__.py b/src/calng/geometries/__init__.py index a13f8291..a2864497 100644 --- a/src/calng/geometries/__init__.py +++ b/src/calng/geometries/__init__.py @@ -6,4 +6,6 @@ from . import ( Lpd1MGeometry, LpdminiGeometry, JungfrauGeometry, + Lpd1MGeometry, + PnccdGeometry, ) diff --git a/src/calng/kernels/common_cpu.pyx b/src/calng/kernels/common_cpu.pyx new file mode 100644 index 00000000..893015f2 --- /dev/null +++ b/src/calng/kernels/common_cpu.pyx @@ -0,0 +1,70 @@ +# cython: boundscheck=False +# cython: cdivision=True +# cython: wrapararound=False + +from cython.parallel import prange +import numpy as np + + +cdef extern from "<algorithm>" namespace "std" nogil: + void nth_element[Iter](Iter first, Iter nth, Iter last) except + + + +def cm_fs( + float[:, :] data, + unsigned char[:, :] mask, + float sigma, + float ratio, +): + # aux will essentially contain the unmasked data + # used in median finding which reorders input, hence copying is needed anyway + cdef float[:, :] aux = np.empty_like(data, order="C") + cdef int num_dark, min_dark, i, j, k + cdef float median + min_dark = <int>(<float>(data.shape[1]) * ratio) + for i in prange(data.shape[0], nogil=True): + num_dark = 0 + for j in range(data.shape[1]): + if mask[i, j] == 0: + aux[i, num_dark] = data[i, j] + num_dark = num_dark + 1 + if num_dark < min_dark: + continue + k = num_dark // 2 + nth_element(&aux[i, 0], (&aux[i, 0]) + k, (&aux[i, 0]) + num_dark) + median = aux[i, k] + if num_dark % 2 == 0: + nth_element(&aux[i, 0], (&aux[i, 0]) + k - 1, (&aux[i, 0]) + num_dark) + median = median / 2 + aux[i, k - 1] / 2 + + for j in range(data.shape[1]): + data[i, j] = data[i, j] - median + + +def cm_ss( + float[:, :] data, + unsigned char[:, :] mask, + float sigma, + float ratio, +): + # transposing to give std::nth_element contiguous input + cdef float[:, ::1] aux = np.empty_like(np.transpose(data), order="C") + cdef int num_dark, min_dark, i, j, k + cdef float median + min_dark = <int>(<float>data.shape[1] * ratio) + for j in prange(data.shape[1], nogil=True): + num_dark = 0 + for i in range(data.shape[0]): + if mask[i, j] == 0: + aux[j, num_dark] = data[i, j] + num_dark = num_dark + 1 + if num_dark < min_dark: + continue + k = num_dark // 2 + nth_element(&aux[j, 0], (&aux[j, 0]) + k, (&aux[j, 0]) + <int>num_dark) + median = aux[j, k] + if num_dark % 2 == 0: + nth_element(&aux[j, 0], (&aux[j, 0]) + k - 1, (&aux[j, 0]) + <int>num_dark) + median = median / 2 + aux[j, k - 1] / 2 + for i in range(data.shape[0]): + data[i, j] = data[i, j] - median diff --git a/src/calng/scenes.py b/src/calng/scenes.py index 3849d682..e48ae88b 100644 --- a/src/calng/scenes.py +++ b/src/calng/scenes.py @@ -865,6 +865,74 @@ class ManualOriginGeometrySettings(VerticalLayout): ) +@titled("Manual geometry settings") +@boxed +class RelativeGeometrySettings(VerticalLayout): + def __init__(self, device_id): + super().__init__(padding=0) + self.children.extend( + [ + HorizontalLayout( + LabelModel(text="Gap", width=4 * BASE_INC, height=BASE_INC), + DoubleLineEditModel( + keys=[f"{device_id}.manualSetting.gap"], + width=4 * BASE_INC, + height=BASE_INC, + ), + ), + HorizontalLayout( + LabelModel( + text="Module offsets", width=6 * BASE_INC, height=BASE_INC + ), + LabelModel(text="x", width=4 * BASE_INC, height=BASE_INC), + LabelModel(text="y", width=4 * BASE_INC, height=BASE_INC), + LabelModel(text="z", width=4 * BASE_INC, height=BASE_INC), + ), + HorizontalLayout( + LabelModel(text="Top", width=6 * BASE_INC, height=BASE_INC), + DoubleLineEditModel( + keys=[f"{device_id}.manualSetting.topOffset.x"], + width=4 * BASE_INC, + height=BASE_INC, + ), + DoubleLineEditModel( + keys=[f"{device_id}.manualSetting.topOffset.y"], + width=4 * BASE_INC, + height=BASE_INC, + ), + DoubleLineEditModel( + keys=[f"{device_id}.manualSetting.topOffset.z"], + width=4 * BASE_INC, + height=BASE_INC, + ), + ), + HorizontalLayout( + LabelModel(text="Bottom", width=6 * BASE_INC, height=BASE_INC), + DoubleLineEditModel( + keys=[f"{device_id}.manualSetting.bottomOffset.x"], + width=4 * BASE_INC, + height=BASE_INC, + ), + DoubleLineEditModel( + keys=[f"{device_id}.manualSetting.bottomOffset.y"], + width=4 * BASE_INC, + height=BASE_INC, + ), + DoubleLineEditModel( + keys=[f"{device_id}.manualSetting.bottomOffset.z"], + width=4 * BASE_INC, + height=BASE_INC, + ), + ), + DisplayCommandModel( + keys=[f"{device_id}.manualSetting.setManual"], + width=6 * BASE_INC, + height=BASE_INC, + ), + ] + ) + + @titled("Manual geometry settings") @boxed class ManualModulesGeometrySettings(VerticalLayout): @@ -1732,7 +1800,7 @@ def quadrant_geometry_overview(device_id, schema): @scene_generator def origin_geometry_overview(device_id, schema): - schema_hash = schema_to_hash(schema) + # TODO: handle loading return VerticalLayout( HorizontalLayout( ManualOriginGeometrySettings(device_id), @@ -1743,6 +1811,18 @@ def origin_geometry_overview(device_id, schema): ) +@scene_generator +def relative_geometry_overview(device_id, schema): + return VerticalLayout( + HorizontalLayout( + RelativeGeometrySettings(device_id), + # GeometryFromFileSettings(device_id, schema_hash), + # TweakCurrentGeometry(device_id), + ), + GeometryPreview(device_id), + ) + + @scene_generator def modules_geometry_overview(device_id, schema): schema_hash = schema_to_hash(schema) diff --git a/src/calng/schemas.py b/src/calng/schemas.py index 2443c9c3..bbe2f9e7 100644 --- a/src/calng/schemas.py +++ b/src/calng/schemas.py @@ -5,6 +5,7 @@ from karabo.bound import ( NDARRAY_ELEMENT, NODE_ELEMENT, STRING_ELEMENT, + UINT64_ELEMENT, VECTOR_STRING_ELEMENT, Schema, ) @@ -256,3 +257,33 @@ def jf_output_schema(use_shmem_handle=True): .commit(), ) return res + + +def pnccd_output_schema(use_shmem_handle=True): + res = Schema() + ( + NODE_ELEMENT(res) + .key("data") + .commit(), + + UINT64_ELEMENT(res) + .key("data.trainId") + .assignmentOptional() + .defaultValue(0) + .commit(), + ) + if use_shmem_handle: + ( + STRING_ELEMENT(res) + .key("data.image") + .assignmentOptional() + .noDefaultValue() + .commit(), + ) + else: + ( + NDARRAY_ELEMENT(res) + .key("data.image") + .commit(), + ) + return res diff --git a/src/tests/test_pnccd_kernels.py b/src/tests/test_pnccd_kernels.py new file mode 100644 index 00000000..18f438fa --- /dev/null +++ b/src/tests/test_pnccd_kernels.py @@ -0,0 +1,50 @@ +import numpy as np +from calng import utils +from calng.corrections.PnccdCorrection import Constants, CorrectionFlags, PnccdCpuRunner + + +def test_common_mode(): + def slow_common_mode( + data, bad_pixel_map, noise_map, cm_min_frac, cm_noise_sigma, cm_ss, cm_fs + ): + masked = np.ma.masked_array( + data=data, + mask=(bad_pixel_map != 0) | (data > noise_map * cm_noise_sigma), + ) + + if cm_fs: + subset_fs = masked.count(axis=1) >= masked.shape[1] * cm_min_frac + data[subset_fs] -= np.ma.median(masked[subset_fs], axis=1, keepdims=True) + + if cm_ss: + subset_ss = masked.count(axis=0) >= masked.shape[0] * cm_min_frac + data[:, subset_ss] -= np.ma.median(masked[:, subset_ss], axis=0) + + data_shape = (1024, 1024) + data = (np.random.random(data_shape) * 1000).astype(np.uint16) + runner = PnccdCpuRunner(1024, 1024, 1, 1) + for noise_level in (0, 5, 100, 1000): + noise = (np.random.random(data_shape) * noise_level).astype(np.float32) + runner.load_constant(Constants.NoiseCCD, noise) + for bad_pixel_percentage in (0, 1, 20, 80, 100): + bad_pixels = ( + np.random.random(data_shape) < (bad_pixel_percentage / 100) + ).astype(np.uint32) + runner.load_constant(Constants.BadPixelsDarkCCD, bad_pixels) + runner.load_data(data.copy()) + runner.correct( + CorrectionFlags.COMMONMODE, + cm_min_frac=0.25, + cm_noise_sigma=10, + cm_row=True, + cm_col=True, + ) + + slow_version = data.astype(np.float32, copy=True) + for q_data, q_noise, q_bad_pixels in zip( + utils.quadrant_views(slow_version), + utils.quadrant_views(noise), + utils.quadrant_views(bad_pixels), + ): + slow_common_mode(q_data, q_bad_pixels, q_noise, 0.25, 10, True, True) + assert np.allclose(runner.processed_data, slow_version, equal_nan=True) -- GitLab