From 366b9071f80ae148ab7f32de4347da9579cb136e Mon Sep 17 00:00:00 2001
From: Oleksii Turkot <oleksii.turkot@xfel.eu>
Date: Tue, 22 Jun 2021 14:49:05 +0100
Subject: [PATCH 1/9] Bring in mask conversion code from EXtra-xwiz
https://git.xfel.eu/gitlab/dataAnalysis/crystfel_spb_workflow/tree/dev/extra_xwiz/mask_converter
Commit ddf45d59f03b7ff2ba7768e70c438baa43e335b5
---
extra_geom/mask.py | 489 +++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 489 insertions(+)
create mode 100644 extra_geom/mask.py
diff --git a/extra_geom/mask.py b/extra_geom/mask.py
new file mode 100644
index 00000000..35e4b47f
--- /dev/null
+++ b/extra_geom/mask.py
@@ -0,0 +1,489 @@
+"""Convert between
+"""
+
+import os
+import re
+import warnings
+
+import h5py
+import numpy as np
+
+
+def delta_method(mask):
+ """
+ Generalized delta-method algorithm to decompose two-dimensional mask
+ (boolean numpy array) into rectangles.
+
+ Based on paper:
+ http://library.utia.cas.cz/separaty/2012/ZOI
+ /suk-rectangular%20decomposition%20of%20binary%20images.pdf
+
+ Args:
+ mask (2D np.array, dtype=bool): input mask.
+
+ Raises:
+ ValueError: mask is not a 2D boolean numpy array.
+
+ Returns:
+ list: list of rectangles that represent masked regionsin in the form:
+ [((x_min, x_max),(y_min, y_max)), ...]
+ """
+
+ # Check input:
+ if (mask.ndim != 2
+ or mask.dtype != bool):
+ raise ValueError(f"Expected input - 2D boolean numpy array.")
+
+ A = np.copy(mask)
+ res = []
+
+ i_y = 0
+ while not np.array_equiv(A, False):
+
+ assert i_y < A.shape[0], (
+ "Decomposition GDM: algorithm did not convert the whole matrix.")
+
+ x_min = x_max = y_min = y_max = -1
+ for i_x in range(A.shape[1]):
+ if x_min < 0 and A[i_y, i_x]:
+ x_min = i_x
+ y_min = i_y
+ if i_x == A.shape[1] - 1:
+ x_max = A.shape[1]
+ elif x_min >= 0 and not A[i_y, i_x]:
+ x_max = i_x
+ elif x_min >= 0 and i_x == A.shape[1] - 1:
+ x_max = A.shape[1]
+
+ if x_max >= 0:
+ sx = slice(x_min, x_max)
+ for i_y2 in range(i_y + 1, A.shape[0]):
+ if not np.array_equal(A[i_y2, sx], A[i_y, sx]):
+ y_max = i_y2
+ break
+ else:
+ y_max = A.shape[0]
+ res.append(((x_min, x_max), (y_min, y_max)))
+ A[y_min:y_max, sx] = False
+ x_min = x_max = y_min = y_max = -1
+
+ i_y += 1
+
+ return res
+
+
+class MaskConverter:
+
+ modes_avail = ('hd52geom', 'geom2hd5')
+ write_modes = ('replace', 'add')
+
+ def __init__(self, hd5file, geofile, run_mode, write_mode,
+ hd5path, hd5entry, detector, data_type, invert):
+ """
+ Construct a mask converter with provided parameters.
+
+ Args:
+ hd5file (string): relative or absolute path to the HD5 file.
+ geofile (string): relative or absolute path to the geometry file.
+ run_mode (string): mode of the mask converter operation.
+ write_mode (string): mode of writing to the existing file.
+ hd5path (string): path to the mask in HD5 file.
+ hd5entry (int): entry number of the mask in HD5 file.
+ detector (string): name of the detector.
+ data_type (string): type of the detector data.
+ invert (bool): invert the mask after reading from of before
+ writing to the HD5 file.
+
+ Raises:
+ KeyError: unexpected run mode.
+ KeyError: unexpected write mode.
+ KeyError: detector name missing in the detector_info keys.
+ KeyError: data type missing in detector_info for specified
+ detector.
+ """
+
+ self._hd5file = hd5file
+ self._hd5path = hd5path
+ self._hd5entry = hd5entry
+ self._geofile = geofile
+ self._run_mode = run_mode
+ self._write_mode = write_mode
+ self._detector = detector
+ self._data_type = data_type
+ self._invert = invert
+
+ if self._run_mode not in self.modes_avail:
+ raise KeyError(
+ f"Unknown MaskConverter run mode: {self._run_mode}.")
+ if self._write_mode not in self.write_modes:
+ raise KeyError(
+ f"Unknown MaskConverter write mode: {self._write_mode}.")
+ if self._detector not in di.detector_info.keys():
+ raise KeyError(
+ f"Unknown detector: {self._detector}.")
+ if self._data_type not in di.detector_info[self._detector].keys():
+ raise KeyError(
+ f"Missing info on '{self._data_type}' data type for "
+ f"detector {self._detector}.")
+
+ self._det_info = di.detector_info[self._detector][self._data_type]
+
+ self._read_mask()
+
+ @property
+ def mask(self):
+ """
+ Provides a copy of the detector mask.
+
+ Returns:
+ np.array: Detector mask as a 2D boolean numpy array.
+ """
+ return np.copy(self.__mask)
+
+ def convert(self):
+ """
+ Convert detector mask and write to the output file.
+ """
+ self._convert_mask()
+ self._write_mask()
+
+ def _read_mask_hd5(self):
+ """
+ Read mask from the HD5 file to boolean numpy array.
+
+ Returns:
+ np.array: Detector mask as a 2D boolean numpy array.
+ """
+
+ res_mask = None
+ mask_shape = self._det_info['shape']
+
+ with h5py.File(self._hd5file, 'r') as f_hd5:
+ hd5mask = f_hd5.get(self._hd5path)
+
+ # Check shape of the HDF5 mask
+ mu.check_hd5mask(hd5mask, mask_shape, self._hd5entry)
+
+ # Remove 'n_data' dimension, convert to np.array
+ if hd5mask.ndim > len(mask_shape):
+ hd5mask = hd5mask[self._hd5entry]
+ else:
+ hd5mask = np.array(hd5mask)
+
+ # Convert mask values to True (masked) and False
+ res_mask = self._det_info['read_mask'](hd5mask)
+
+ if self._invert:
+ res_mask = np.logical_not(res_mask)
+
+ return res_mask
+
+ def _read_mask_geo(self):
+ """
+ Read mask from the geometry file to the dictionary of rectangles
+ (in the same format as in the geometry file).
+
+ Raises:
+ ValueError: mask panel description not expected for the
+ specified data type.
+ ValueError: mask panel description does not suit specified
+ detector and data type.
+
+ Returns:
+ dict: dictionary of rectangles (in the same format as in the
+ geometry file) representing masked regions.
+ """
+
+ res_dict = {}
+ bad_dict = {}
+
+ with open(self._geofile, 'r') as f_geo:
+ for line in f_geo:
+ # To ignore comments:
+ line = line.partition(';')[0].rstrip()
+
+ m_obj = re.search(r"bad_(.+)/(\S+)\s*=\s*(\S+)", line)
+ if m_obj is not None:
+ area, var, val = m_obj.groups()
+ if area not in bad_dict.keys():
+ bad_dict[area] = {
+ 'min_fs': -1,
+ 'max_fs': -1,
+ 'min_ss': -1,
+ 'max_ss': -1,
+ 'panel': 'all'
+ }
+ if var == 'panel':
+ bad_dict[area][var] = val
+
+ # Check whether panel value is expected
+ if len(self._det_info['shape']) != 3:
+ raise ValueError(
+ f"Panel description ({val}) not expected "
+ f"for {self._data_type} data.")
+
+ # Check <val> to be suitable description
+ # of a panel and asic
+ m_obj = re.match(r"(p\d+)(a\d+)", val)
+ if (m_obj is None
+ or m_obj.group(1) not in
+ self._det_info['panel_names']
+ or m_obj.group(2) not in
+ self._det_info['asic_names']):
+ raise ValueError(
+ f"Not suitable panel description: {val}.")
+ elif var in bad_dict[area].keys():
+ bad_dict[area][var] = int(val)
+ else:
+ warnings.warn(
+ f"Geometry file - unsupported mask variable: "
+ f"{var} in bad_{area}.")
+
+ # Check if rectangle information is complete
+ i = 0
+ for area in bad_dict.keys():
+ if all([
+ bad_dict[area]['min_fs'] >= 0,
+ bad_dict[area]['max_fs'] >= 0,
+ bad_dict[area]['min_ss'] >= 0,
+ bad_dict[area]['max_ss'] >= 0
+ ]):
+ res_dict[i] = bad_dict[area]
+ i += 1
+ else:
+ warnings.warn(
+ f"Geometry file - incomplete information for bad_{area}.")
+
+ return res_dict
+
+ def _read_mask(self):
+ """
+ Read mask from the HD5 or geometry file, depending on mode.
+ """
+
+ self.__mask = np.zeros(self._det_info['shape'], dtype=bool)
+ self.__rect = {}
+
+ if self._run_mode == "hd52geom":
+ self.__mask = self._read_mask_hd5()
+
+ # Reduce mask in case of write option 'add'
+ if self._write_mode == 'add' and os.path.exists(self._geofile):
+ self.__rect = self._read_mask_geo()
+ reduce_mask = self._convert_rectd2nparr()
+ self.__mask = np.logical_and(self.__mask,
+ np.logical_not(reduce_mask))
+ self.__rect = {}
+
+ elif self._run_mode == "geom2hd5":
+ self.__rect = self._read_mask_geo()
+
+ # Read also mask from HD5 in case of write option 'add'
+ if (self._write_mode == 'add'
+ and os.path.exists(self._hd5file)):
+ self.__mask = self._read_mask_hd5()
+
+ def _convert_nparr2rectd(self):
+ """
+ Convert mask from the 2D boolean numpy array to the dictionary
+ of rectangles (in the same format as in the geometry file).
+
+ Returns:
+ dict: Dictionary of rectangles (in the same format as in the
+ geometry file) representing masked regions.
+ """
+
+ res_dict = {}
+
+ # First check for regions to be excluded in all panels:
+ if self.__mask.ndim == 3:
+ panel_all = self.mask[0]
+ for i in range(1, self.__mask.shape[0]):
+ panel_all = np.logical_and(panel_all, self.__mask[i])
+ else:
+ panel_all = self.mask
+
+ is_panel_all_empty = np.array_equiv(panel_all, False)
+ if not is_panel_all_empty:
+ res_dict.update(mu.rect2dict(delta_method(panel_all), 'all'))
+
+ if self.__mask.ndim == 3:
+ panels = self._det_info['panel_names']
+ asics = self._det_info['asic_names']
+ asic_range = self._det_info['asic_range']
+
+ # Loop over all panels:
+ for i in range(len(panels)):
+ if is_panel_all_empty:
+ panel_i = self.__mask[i]
+ else:
+ panel_i = np.copy(self.__mask[i])
+ panel_i = np.logical_and(panel_i,
+ np.logical_not(panel_all))
+
+ # Loop over all asics in the panel:
+ for j in range(len(asics)):
+ asic_j = np.zeros(panel_i.shape, dtype=bool)
+ asic_j[asic_range[i][j]] = True
+ panel_i_asic_j = np.logical_and(panel_i, asic_j)
+ res_dict.update(
+ mu.rect2dict(delta_method(panel_i_asic_j),
+ f"{panels[i]}{asics[j]}"))
+
+ return res_dict
+
+ def _convert_rectd2nparr(self):
+ """
+ Convert mask from the dictionary of rectangles (same format as
+ in the geometry file) to the 2D boolean numpy array.
+
+ Returns:
+ np.array: Detector mask as a 2D boolean numpy array.
+ """
+
+ shape = self._det_info['shape']
+ res_mask = np.zeros(shape, dtype=bool)
+
+ for area in self.__rect.keys():
+ slice_ss = slice(self.__rect[area]['min_ss'],
+ self.__rect[area]['max_ss'] + 1)
+ slice_fs = slice(self.__rect[area]['min_fs'],
+ self.__rect[area]['max_fs'] + 1)
+ if self.__rect[area]['panel'] == 'all':
+ if len(shape) == 3:
+ res_mask[:, slice_ss, slice_fs] = True
+ else:
+ res_mask[slice_ss, slice_fs] = True
+ else:
+ assert len(shape) == 3, (
+ "Convert rectd2nparr: mask has to be dimensions 3 to "
+ "apply rectangles per panel.")
+
+ panel, asic = re.match(
+ r"(p\d+)(a\d+)", self.__rect[area]['panel']).groups()
+ panel_n = self._det_info['panel_names'].index(panel)
+ asic_n = self._det_info['asic_names'].index(asic)
+ asic_range = self._det_info['asic_range'][panel_n][asic_n]
+
+ slice_ss_in_asic = slice(
+ max(slice_ss.start, asic_range[0].start),
+ min(slice_ss.stop, asic_range[0].stop))
+ slice_fs_in_asic = slice(
+ max(slice_fs.start, asic_range[1].start),
+ min(slice_fs.stop, asic_range[1].stop))
+
+ res_mask[panel_n, slice_ss_in_asic, slice_fs_in_asic] = True
+
+ return res_mask
+
+ def _convert_mask(self):
+ """
+ Convert the mask, depending on mode.
+ """
+ if self._run_mode == "hd52geom":
+ self.__rect = self._convert_nparr2rectd()
+ elif self._run_mode == "geom2hd5":
+ rect_mask = self._convert_rectd2nparr()
+ self.__mask = np.logical_or(self.__mask, rect_mask)
+
+ def _write_mask_hd5(self):
+ """
+ Write converted mask to the HD5 file.
+ """
+
+ mask_shape = self.__mask.shape
+ mask_tmp = self.__mask
+ if self._invert:
+ mask_tmp = np.logical_not(mask_tmp)
+ mask_to_write = self._det_info['write_mask'](mask_tmp)
+
+ with h5py.File(self._hd5file, 'a') as f_hd5:
+ if self._hd5path in f_hd5:
+ hd5mask = f_hd5[self._hd5path]
+
+ # Check shape of the existing HDF5 mask
+ mu.check_hd5mask(hd5mask, mask_shape, self._hd5entry)
+
+ if hd5mask.ndim > len(mask_shape):
+ hd5mask[self._hd5entry] = mask_to_write
+ else:
+ hd5mask[...] = mask_to_write
+ else:
+ f_hd5.create_dataset(self._hd5path, data=mask_to_write)
+
+ def _write_mask_geo(self):
+ """
+ Write converted mask to the geometry file.
+ """
+
+ text_before = text_after = []
+ n_area_start = 0
+
+ # Store and process content of the existing geometry file
+ if os.path.exists(self._geofile):
+ with open(self._geofile, 'r') as f_geo:
+ contents = f_geo.readlines()
+
+ idx_write = len(contents)
+ for i, line in enumerate(contents):
+
+ # Search for the mask information
+ if "bad_" in line:
+ idx_write = i + 1
+
+ # Comment existing mask for the 'replace' mode
+ if all([
+ self._write_mode == 'replace',
+ "bad_" in line.partition(';')[0]
+ ]):
+ contents[i] = "; " + line
+
+ # Check existing numbered 'bad_area's
+ sa_obj = re.search(r"bad_area(\d+)/", line)
+ if sa_obj is not None:
+ sa_num = int(sa_obj.group(1))
+ n_area_start = max(n_area_start, sa_num + 1)
+
+ # Search for the 'rigid_group'
+ # (to put mask before it in case no 'bad_' area found)
+ if all([
+ idx_write == len(contents),
+ "rigid_group" in line
+ ]):
+ idx_write = i
+
+ text_before = contents[:idx_write]
+ text_after = contents[idx_write:]
+
+ # Format mask as a list of text lines
+ text_mask = []
+ if (text_before
+ and text_before[-1].strip() != ""):
+ text_mask.append("\n")
+
+ for area in self.__rect:
+ for dim in ['min_fs', 'max_fs', 'min_ss', 'max_ss']:
+ text_mask.append(
+ f"bad_area{n_area_start + area}/{dim} = "
+ f"{self.__rect[area][dim]}\n")
+ if self.__rect[area]['panel'] != 'all':
+ text_mask.append(
+ f"bad_area{n_area_start + area}/panel = "
+ f"{self.__rect[area]['panel']}\n")
+ text_mask.append("\n")
+
+ text_write = "".join(text_before + text_mask + text_after)
+
+ with open(self._geofile, 'w') as f_geo:
+ f_geo.write(text_write)
+
+ def _write_mask(self):
+ """
+ Write converted mask to the HD5 or geometry file, depending on mode.
+ """
+ # Write the mask depending on mode.
+ if self._run_mode == "hd52geom":
+ self._write_mask_geo()
+ elif self._run_mode == "geom2hd5":
+ self._write_mask_hd5()
+
From fc4c37bd36297e39831fec26e5f677d0a5a64c9a Mon Sep 17 00:00:00 2001
From: Thomas Kluyver <thomas@kluyver.me.uk>
Date: Wed, 23 Jun 2021 14:28:24 +0100
Subject: [PATCH 2/9] Refactoring mask conversion code
---
extra_geom/crystfel_fmt.py | 43 +++-
extra_geom/mask.py | 508 +++++++++++++------------------------
2 files changed, 214 insertions(+), 337 deletions(-)
diff --git a/extra_geom/crystfel_fmt.py b/extra_geom/crystfel_fmt.py
index 540875d8..cbfb8af4 100644
--- a/extra_geom/crystfel_fmt.py
+++ b/extra_geom/crystfel_fmt.py
@@ -1,4 +1,4 @@
-"""Write geometry in CrystFEL format.
+"""Write & process geometry in CrystFEL format.
"""
from itertools import product
@@ -190,3 +190,44 @@ def get_rigid_groups(geom, nquads=4):
]
return '\n'.join(lines)
+
+
+def panel_modno(panel_info, pname):
+ """Find the module number (counting from 0) of the given panel
+
+ Returns None if the geometry describes a 2D input array rather than a 3D
+ array with a stack of modules.
+ """
+ dims = panel_info['dim_structure']
+ ix_dims = [i for i in dims if isinstance(i, int)]
+ if len(ix_dims) > 1:
+ raise ValueError(f"Too many index dimensions for {pname}: {dims}")
+
+ return ix_dims[0] if ix_dims else None
+
+
+def data_shape(panels_dict):
+ """Make a 2- or 3-tuple representing the expected data shape
+
+ panels_dict is the structure given by cfelpytuils::
+
+ load_crystfel_geometry(f)['panels']
+ """
+ nmodules = slow_scan_px = fast_scan_px = 0
+
+ for pname, info in panels_dict.items():
+ modno = panel_modno(info, pname)
+ if modno is not None:
+ nmodules = max(nmodules, modno + 1)
+ slow_scan_px = max(slow_scan_px, info['max_ss'] + 1)
+ fast_scan_px = max(fast_scan_px, info['max_fs'] + 1)
+
+ if nmodules == 0:
+ shape = (slow_scan_px, fast_scan_px)
+ else:
+ shape = (nmodules, slow_scan_px, fast_scan_px)
+
+ if min(shape) <= 0:
+ raise ValueError("Could not find detector data shape from .geom file")
+
+ return shape
diff --git a/extra_geom/mask.py b/extra_geom/mask.py
index 35e4b47f..e01297c4 100644
--- a/extra_geom/mask.py
+++ b/extra_geom/mask.py
@@ -1,13 +1,15 @@
-"""Convert between
+"""Convert between masks as arrays and as a set of rectangular regions
"""
import os
import re
import warnings
-import h5py
+from cfelpyutils.crystfel_utils import load_crystfel_geometry
import numpy as np
+from . import crystfel_fmt
+
def delta_method(mask):
"""
@@ -71,357 +73,193 @@ def delta_method(mask):
return res
+class RegionRect:
+ """One masked region within a 3D array of data (modules stacked)
-class MaskConverter:
+ This can apply to a single module (modno=5) or all modules (modno=None).
+ """
+ def __init__(self, modno, start_ss, stop_ss, start_fs, stop_fs):
+ self.modno = modno # None -> applies to all modules
+ self.start_ss = start_ss
+ self.stop_ss = stop_ss
+ self.start_fs = start_fs
+ self.stop_fs = stop_fs
- modes_avail = ('hd52geom', 'geom2hd5')
- write_modes = ('replace', 'add')
+ def _tuple(self):
+ return self.modno, self.start_ss, self.stop_ss, self.start_fs, self.stop_fs
- def __init__(self, hd5file, geofile, run_mode, write_mode,
- hd5path, hd5entry, detector, data_type, invert):
- """
- Construct a mask converter with provided parameters.
-
- Args:
- hd5file (string): relative or absolute path to the HD5 file.
- geofile (string): relative or absolute path to the geometry file.
- run_mode (string): mode of the mask converter operation.
- write_mode (string): mode of writing to the existing file.
- hd5path (string): path to the mask in HD5 file.
- hd5entry (int): entry number of the mask in HD5 file.
- detector (string): name of the detector.
- data_type (string): type of the detector data.
- invert (bool): invert the mask after reading from of before
- writing to the HD5 file.
-
- Raises:
- KeyError: unexpected run mode.
- KeyError: unexpected write mode.
- KeyError: detector name missing in the detector_info keys.
- KeyError: data type missing in detector_info for specified
- detector.
- """
+ def __hash__(self):
+ return hash(self._tuple())
- self._hd5file = hd5file
- self._hd5path = hd5path
- self._hd5entry = hd5entry
- self._geofile = geofile
- self._run_mode = run_mode
- self._write_mode = write_mode
- self._detector = detector
- self._data_type = data_type
- self._invert = invert
-
- if self._run_mode not in self.modes_avail:
- raise KeyError(
- f"Unknown MaskConverter run mode: {self._run_mode}.")
- if self._write_mode not in self.write_modes:
- raise KeyError(
- f"Unknown MaskConverter write mode: {self._write_mode}.")
- if self._detector not in di.detector_info.keys():
- raise KeyError(
- f"Unknown detector: {self._detector}.")
- if self._data_type not in di.detector_info[self._detector].keys():
- raise KeyError(
- f"Missing info on '{self._data_type}' data type for "
- f"detector {self._detector}.")
-
- self._det_info = di.detector_info[self._detector][self._data_type]
-
- self._read_mask()
+ def __eq__(self, other):
+ return isinstance(other, RegionRect) and self._tuple() == other._tuple()
@property
- def mask(self):
- """
- Provides a copy of the detector mask.
+ def array_slice(self):
+ """Get a tuple to use for slicing an array"""
+ mod = np.s_[:] if (self.modno is None) else self.modno
+ return np.s_[mod, self.start_ss:self.stop_ss, self.start_fs:self.stop_fs]
- Returns:
- np.array: Detector mask as a 2D boolean numpy array.
- """
- return np.copy(self.__mask)
+ def intersection(self, other: 'RegionRect'):
+ """Find the intersection of this and another RegionRect
- def convert(self):
+ Returns (overlap, intersection). The first value is True if the two
+ regions overlap, and the second is a RegionRect for their intersection.
"""
- Convert detector mask and write to the output file.
+ modno = self.modno if (self.modno is not None) else other.modno
+ if other.modno in (modno, None):
+ # On the same module, or one of self/other is for all modules
+ r = RegionRect(
+ modno,
+ max(self.start_ss, other.start_ss),
+ min(self.stop_ss, other.stop_ss),
+ max(self.start_fs, other.start_fs),
+ min(self.stop_fs, other.stop_fs),
+ )
+ overlap = (r.start_ss < r.stop_ss) and (r.start_fs < r.stop_fs)
+ return overlap, r
+
+ return False, RegionRect(0, 0, 0, 0, 0)
+
+
+class MaskRegions:
+ """A set of rectangular regions masked in a stack of module images"""
+ def __init__(self, shape, regions=()):
+ assert len(shape) == 3
+ self.shape = shape
+ self.regions = regions # [RegionRect]
+
+ @classmethod
+ def from_crystfel_geom(cls, filename):
+ """Read mask from a CrystFEL format geometry file.
"""
- self._convert_mask()
- self._write_mask()
-
- def _read_mask_hd5(self):
- """
- Read mask from the HD5 file to boolean numpy array.
-
- Returns:
- np.array: Detector mask as a 2D boolean numpy array.
- """
-
- res_mask = None
- mask_shape = self._det_info['shape']
-
- with h5py.File(self._hd5file, 'r') as f_hd5:
- hd5mask = f_hd5.get(self._hd5path)
-
- # Check shape of the HDF5 mask
- mu.check_hd5mask(hd5mask, mask_shape, self._hd5entry)
-
- # Remove 'n_data' dimension, convert to np.array
- if hd5mask.ndim > len(mask_shape):
- hd5mask = hd5mask[self._hd5entry]
+ geom_dict = load_crystfel_geometry(filename)
+ arr_shape = crystfel_fmt.data_shape(geom_dict['panels'])
+ if len(arr_shape) == 2:
+ arr_shape = (1,) + arr_shape
+
+ regions = []
+ for name, info in geom_dict['bad'].items():
+ if not info['is_fsss']:
+ # Regions can also be defined with x & y coordinates, but these
+ # are not handled yet.
+ continue
+
+ min_ss, max_ss = info['min_ss'], info['max_ss']
+ min_fs, max_fs = info['min_fs'], info['max_fs']
+ if (max_ss < min_ss) or (max_fs < min_fs):
+ warnings.warn(
+ f"Geometry file {filename} - incomplete information for {name}."
+ )
+ if 'panel' in info:
+ panel_info = geom_dict['panels'][info['panel']]
+ modno = crystfel_fmt.panel_modno(panel_info, info['panel'])
+ modno = modno or 0 # None -> 0
+
+ # Clip the region to the panel limits
+ min_ss = max(min_ss, panel_info['min_ss'])
+ max_ss = min(max_ss, panel_info['max_ss'])
+ min_fs = max(min_fs, panel_info['min_fs'])
+ max_fs = min(max_fs, panel_info['max_fs'])
else:
- hd5mask = np.array(hd5mask)
+ modno = None # Applies to all panels
- # Convert mask values to True (masked) and False
- res_mask = self._det_info['read_mask'](hd5mask)
+ regions.append(RegionRect(
+ modno, min_ss, max_ss + 1, min_fs, max_fs + 1,
+ ))
- if self._invert:
- res_mask = np.logical_not(res_mask)
+ return cls(arr_shape, regions)
- return res_mask
+ @classmethod
+ def from_mask_array(cls, arr):
+ """Convert mask from a 2D or 3D boolean numpy array to rectangular regions.
- def _read_mask_geo(self):
- """
- Read mask from the geometry file to the dictionary of rectangles
- (in the same format as in the geometry file).
-
- Raises:
- ValueError: mask panel description not expected for the
- specified data type.
- ValueError: mask panel description does not suit specified
- detector and data type.
-
- Returns:
- dict: dictionary of rectangles (in the same format as in the
- geometry file) representing masked regions.
+ The masked/selected regions will be where the array contains True.
"""
+ if arr.ndim not in {2, 3}:
+ raise TypeError(f"Array must be 2D or 3D (got {arr.ndim} dims)")
- res_dict = {}
- bad_dict = {}
-
- with open(self._geofile, 'r') as f_geo:
- for line in f_geo:
- # To ignore comments:
- line = line.partition(';')[0].rstrip()
-
- m_obj = re.search(r"bad_(.+)/(\S+)\s*=\s*(\S+)", line)
- if m_obj is not None:
- area, var, val = m_obj.groups()
- if area not in bad_dict.keys():
- bad_dict[area] = {
- 'min_fs': -1,
- 'max_fs': -1,
- 'min_ss': -1,
- 'max_ss': -1,
- 'panel': 'all'
- }
- if var == 'panel':
- bad_dict[area][var] = val
-
- # Check whether panel value is expected
- if len(self._det_info['shape']) != 3:
- raise ValueError(
- f"Panel description ({val}) not expected "
- f"for {self._data_type} data.")
-
- # Check <val> to be suitable description
- # of a panel and asic
- m_obj = re.match(r"(p\d+)(a\d+)", val)
- if (m_obj is None
- or m_obj.group(1) not in
- self._det_info['panel_names']
- or m_obj.group(2) not in
- self._det_info['asic_names']):
- raise ValueError(
- f"Not suitable panel description: {val}.")
- elif var in bad_dict[area].keys():
- bad_dict[area][var] = int(val)
- else:
- warnings.warn(
- f"Geometry file - unsupported mask variable: "
- f"{var} in bad_{area}.")
-
- # Check if rectangle information is complete
- i = 0
- for area in bad_dict.keys():
- if all([
- bad_dict[area]['min_fs'] >= 0,
- bad_dict[area]['max_fs'] >= 0,
- bad_dict[area]['min_ss'] >= 0,
- bad_dict[area]['max_ss'] >= 0
- ]):
- res_dict[i] = bad_dict[area]
- i += 1
- else:
- warnings.warn(
- f"Geometry file - incomplete information for bad_{area}.")
+ def find_regions(arr_2d, modno):
+ return [RegionRect(
+ modno, min_ss, max_ss + 1, min_fs, max_fs + 1
+ ) for ((min_fs, max_fs), (min_ss, max_ss)) in delta_method(arr_2d)]
- return res_dict
+ if arr.ndim == 2:
+ return cls((1,) + arr.shape, find_regions(arr, modno=np.s_[:]))
- def _read_mask(self):
- """
- Read mask from the HD5 or geometry file, depending on mode.
- """
-
- self.__mask = np.zeros(self._det_info['shape'], dtype=bool)
- self.__rect = {}
+ # 3D array (modno, slow_scan, fast_scan)
- if self._run_mode == "hd52geom":
- self.__mask = self._read_mask_hd5()
+ # First check for regions to be excluded in all panels:
+ panel_all = np.logical_and.reduce(arr, axis=0)
+ regions = find_regions(panel_all, modno=np.s_[:])
+ is_panel_all_empty = not np.any(panel_all)
- # Reduce mask in case of write option 'add'
- if self._write_mode == 'add' and os.path.exists(self._geofile):
- self.__rect = self._read_mask_geo()
- reduce_mask = self._convert_rectd2nparr()
- self.__mask = np.logical_and(self.__mask,
- np.logical_not(reduce_mask))
- self.__rect = {}
+ # Loop over all panels:
+ for i, panel_arr in enumerate(arr):
+ if not is_panel_all_empty:
+ panel_arr = np.copy(panel_arr) & ~panel_all
- elif self._run_mode == "geom2hd5":
- self.__rect = self._read_mask_geo()
+ regions.extend(find_regions(panel_arr, modno=i))
- # Read also mask from HD5 in case of write option 'add'
- if (self._write_mode == 'add'
- and os.path.exists(self._hd5file)):
- self.__mask = self._read_mask_hd5()
+ return cls(arr.shape, regions)
- def _convert_nparr2rectd(self):
- """
- Convert mask from the 2D boolean numpy array to the dictionary
- of rectangles (in the same format as in the geometry file).
+ def to_mask_array(self):
+ """Convert the mask rectangles to a 3D boolean numpy array.
- Returns:
- dict: Dictionary of rectangles (in the same format as in the
- geometry file) representing masked regions.
+ The selected regions will be True in the array.
"""
+ res_mask = np.zeros(self.shape, dtype=np.bool_)
- res_dict = {}
-
- # First check for regions to be excluded in all panels:
- if self.__mask.ndim == 3:
- panel_all = self.mask[0]
- for i in range(1, self.__mask.shape[0]):
- panel_all = np.logical_and(panel_all, self.__mask[i])
- else:
- panel_all = self.mask
-
- is_panel_all_empty = np.array_equiv(panel_all, False)
- if not is_panel_all_empty:
- res_dict.update(mu.rect2dict(delta_method(panel_all), 'all'))
-
- if self.__mask.ndim == 3:
- panels = self._det_info['panel_names']
- asics = self._det_info['asic_names']
- asic_range = self._det_info['asic_range']
-
- # Loop over all panels:
- for i in range(len(panels)):
- if is_panel_all_empty:
- panel_i = self.__mask[i]
- else:
- panel_i = np.copy(self.__mask[i])
- panel_i = np.logical_and(panel_i,
- np.logical_not(panel_all))
-
- # Loop over all asics in the panel:
- for j in range(len(asics)):
- asic_j = np.zeros(panel_i.shape, dtype=bool)
- asic_j[asic_range[i][j]] = True
- panel_i_asic_j = np.logical_and(panel_i, asic_j)
- res_dict.update(
- mu.rect2dict(delta_method(panel_i_asic_j),
- f"{panels[i]}{asics[j]}"))
-
- return res_dict
-
- def _convert_rectd2nparr(self):
- """
- Convert mask from the dictionary of rectangles (same format as
- in the geometry file) to the 2D boolean numpy array.
-
- Returns:
- np.array: Detector mask as a 2D boolean numpy array.
- """
-
- shape = self._det_info['shape']
- res_mask = np.zeros(shape, dtype=bool)
-
- for area in self.__rect.keys():
- slice_ss = slice(self.__rect[area]['min_ss'],
- self.__rect[area]['max_ss'] + 1)
- slice_fs = slice(self.__rect[area]['min_fs'],
- self.__rect[area]['max_fs'] + 1)
- if self.__rect[area]['panel'] == 'all':
- if len(shape) == 3:
- res_mask[:, slice_ss, slice_fs] = True
- else:
- res_mask[slice_ss, slice_fs] = True
- else:
- assert len(shape) == 3, (
- "Convert rectd2nparr: mask has to be dimensions 3 to "
- "apply rectangles per panel.")
-
- panel, asic = re.match(
- r"(p\d+)(a\d+)", self.__rect[area]['panel']).groups()
- panel_n = self._det_info['panel_names'].index(panel)
- asic_n = self._det_info['asic_names'].index(asic)
- asic_range = self._det_info['asic_range'][panel_n][asic_n]
-
- slice_ss_in_asic = slice(
- max(slice_ss.start, asic_range[0].start),
- min(slice_ss.stop, asic_range[0].stop))
- slice_fs_in_asic = slice(
- max(slice_fs.start, asic_range[1].start),
- min(slice_fs.stop, asic_range[1].stop))
-
- res_mask[panel_n, slice_ss_in_asic, slice_fs_in_asic] = True
+ for region in self.regions:
+ res_mask[region.array_slice] = True
return res_mask
- def _convert_mask(self):
- """
- Convert the mask, depending on mode.
- """
- if self._run_mode == "hd52geom":
- self.__rect = self._convert_nparr2rectd()
- elif self._run_mode == "geom2hd5":
- rect_mask = self._convert_rectd2nparr()
- self.__mask = np.logical_or(self.__mask, rect_mask)
+ def make_crystfel_bad_regions(self, panels_dict):
+ modno_to_panels = {}
+ for pname, pinfo in panels_dict:
+ modno = crystfel_fmt.panel_modno(pinfo, pname)
+ modno_to_panels.setdefault(modno, []).append((pname, pinfo))
- def _write_mask_hd5(self):
- """
- Write converted mask to the HD5 file.
- """
-
- mask_shape = self.__mask.shape
- mask_tmp = self.__mask
- if self._invert:
- mask_tmp = np.logical_not(mask_tmp)
- mask_to_write = self._det_info['write_mask'](mask_tmp)
-
- with h5py.File(self._hd5file, 'a') as f_hd5:
- if self._hd5path in f_hd5:
- hd5mask = f_hd5[self._hd5path]
+ def to_dict(region: RegionRect):
+ return {
+ 'min_ss': region.start_ss, 'max_ss': region.stop_ss - 1,
+ 'min_fs': region.start_fs, 'max_fs': region.stop_fs - 1,
+ }
- # Check shape of the existing HDF5 mask
- mu.check_hd5mask(hd5mask, mask_shape, self._hd5entry)
+ res = []
+ for mask_region in self.regions:
+ if mask_region.modno is None:
+ res.append(to_dict(mask_region)) # Mask for all panels
- if hd5mask.ndim > len(mask_shape):
- hd5mask[self._hd5entry] = mask_to_write
- else:
- hd5mask[...] = mask_to_write
else:
- f_hd5.create_dataset(self._hd5path, data=mask_to_write)
-
- def _write_mask_geo(self):
- """
- Write converted mask to the geometry file.
+ for pname, pinfo in modno_to_panels[mask_region.modno]:
+ panel_region = RegionRect(
+ mask_region.modno,
+ pinfo['min_ss'], pinfo['max_ss'] + 1,
+ pinfo['min_fs'], pinfo['max_fs'] + 1,
+ )
+ overlap, matching = mask_region.intersection(panel_region)
+ if overlap:
+ d = to_dict(matching)
+ d['panel'] = pname
+ res.append(d)
+
+ return res
+
+ def write_crystfel_geom(self, filename, write_mode):
+ """Write the mask regions to a CrystFEL format geometry file.
+
+ *filename* needs to exist already; it will be read to identify the
+ panels which the mask applies to.
"""
text_before = text_after = []
n_area_start = 0
# Store and process content of the existing geometry file
- if os.path.exists(self._geofile):
- with open(self._geofile, 'r') as f_geo:
+ if os.path.exists(filename):
+ with open(filename, 'r') as f_geo:
contents = f_geo.readlines()
idx_write = len(contents)
@@ -433,7 +271,7 @@ def _write_mask_geo(self):
# Comment existing mask for the 'replace' mode
if all([
- self._write_mode == 'replace',
+ write_mode == 'replace',
"bad_" in line.partition(';')[0]
]):
contents[i] = "; " + line
@@ -457,33 +295,31 @@ def _write_mask_geo(self):
# Format mask as a list of text lines
text_mask = []
- if (text_before
- and text_before[-1].strip() != ""):
+ if text_before and text_before[-1].strip() != "":
text_mask.append("\n")
- for area in self.__rect:
- for dim in ['min_fs', 'max_fs', 'min_ss', 'max_ss']:
- text_mask.append(
- f"bad_area{n_area_start + area}/{dim} = "
- f"{self.__rect[area][dim]}\n")
- if self.__rect[area]['panel'] != 'all':
- text_mask.append(
- f"bad_area{n_area_start + area}/panel = "
- f"{self.__rect[area]['panel']}\n")
- text_mask.append("\n")
+ geom_dict = load_crystfel_geometry(filename)
+
+ new_bad_regions = self.make_crystfel_bad_regions(geom_dict['panels'])
+ for i, bad_d in enumerate(new_bad_regions, start=n_area_start):
+ text_mask.extend([
+ f'bad_area{i}/{k} = {v}' for (k, v) in bad_d.items()
+ ] + ['\n'])
text_write = "".join(text_before + text_mask + text_after)
- with open(self._geofile, 'w') as f_geo:
+ with open(filename, 'w') as f_geo:
f_geo.write(text_write)
- def _write_mask(self):
- """
- Write converted mask to the HD5 or geometry file, depending on mode.
- """
- # Write the mask depending on mode.
- if self._run_mode == "hd52geom":
- self._write_mask_geo()
- elif self._run_mode == "geom2hd5":
- self._write_mask_hd5()
+ def __eq__(self, other):
+ if not (isinstance(other, MaskRegions) and (self.shape == other.shape)):
+ return False
+
+ if set(self.regions) == set(other.regions):
+ return True
+
+ # Two equivalent masks could have regions described in different ways.
+ # Converting them both to arrays is the easiest way to normalise them,
+ # though probably not the most efficient.
+ return np.array_equal(self.to_mask_array(), other.to_mask_array())
From fbeccf7e26ee5b96df44c9653b2f403e7c14cc4d Mon Sep 17 00:00:00 2001
From: Thomas Kluyver <thomas@kluyver.me.uk>
Date: Wed, 23 Jun 2021 17:36:04 +0100
Subject: [PATCH 3/9] Fix some leftover uses of np.s_[:]
---
extra_geom/mask.py | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/extra_geom/mask.py b/extra_geom/mask.py
index e01297c4..ff9b4931 100644
--- a/extra_geom/mask.py
+++ b/extra_geom/mask.py
@@ -185,13 +185,13 @@ def find_regions(arr_2d, modno):
) for ((min_fs, max_fs), (min_ss, max_ss)) in delta_method(arr_2d)]
if arr.ndim == 2:
- return cls((1,) + arr.shape, find_regions(arr, modno=np.s_[:]))
+ return cls((1,) + arr.shape, find_regions(arr, modno=None))
# 3D array (modno, slow_scan, fast_scan)
# First check for regions to be excluded in all panels:
panel_all = np.logical_and.reduce(arr, axis=0)
- regions = find_regions(panel_all, modno=np.s_[:])
+ regions = find_regions(panel_all, modno=None)
is_panel_all_empty = not np.any(panel_all)
# Loop over all panels:
From 4dc579a78d7b7c921788027268e87edab08c94c7 Mon Sep 17 00:00:00 2001
From: Oleksii Turkot <oleksii.turkot@xfel.eu>
Date: Thu, 24 Jun 2021 11:22:49 +0100
Subject: [PATCH 4/9] Add test for delta_method()
---
extra_geom/tests/test_mask.py | 22 ++++++++++++++++++++++
1 file changed, 22 insertions(+)
create mode 100644 extra_geom/tests/test_mask.py
diff --git a/extra_geom/tests/test_mask.py b/extra_geom/tests/test_mask.py
new file mode 100644
index 00000000..20317541
--- /dev/null
+++ b/extra_geom/tests/test_mask.py
@@ -0,0 +1,22 @@
+import numpy as np
+
+from extra_geom import mask
+
+def test_delta_method():
+ A = np.zeros((3, 4), dtype=bool)
+ A[1:2, 1:3] = True
+ B = np.copy(A)
+ assert mask.delta_method(A) == [((1, 3), (1, 2))], "Simple rectangle test."
+ assert np.array_equal(A, B), (
+ "Make sure function does not change the original matrix.")
+
+ C = np.zeros((5, 5), dtype=bool)
+ C[0:3, 1:3] = C[1:3, 0:4] = C[2:4, 2:5] = C[4:5, 1:4] = True
+ exp_res_C = [((1, 3), (0, 3)), ((0, 1), (1, 3)), ((3, 4), (1, 5)),
+ ((4, 5), (2, 4)), ((2, 3), (3, 5)), ((1, 2), (4, 5))]
+ assert mask.delta_method(C) == exp_res_C, "Complex shape test."
+
+ D = np.zeros((5, 5), dtype=bool)
+ D[0:1, 0:1] = D[2:3, 2:3] = True
+ exp_res_D = [((0, 1), (0, 1)), ((2, 3), (2, 3))]
+ assert mask.delta_method(D) == exp_res_D, "Test on separate pixels."
From efea9450a316d41bad9f675cecb848e3f0cb8974 Mon Sep 17 00:00:00 2001
From: turkot <oleksii.turkot@xfel.eu>
Date: Fri, 25 Jun 2021 13:21:39 +0200
Subject: [PATCH 5/9] Fix few small bugs with writing a mask to the geometry
file and using delta_method.
---
extra_geom/mask.py | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/extra_geom/mask.py b/extra_geom/mask.py
index ff9b4931..08c4e7d2 100644
--- a/extra_geom/mask.py
+++ b/extra_geom/mask.py
@@ -181,7 +181,7 @@ def from_mask_array(cls, arr):
def find_regions(arr_2d, modno):
return [RegionRect(
- modno, min_ss, max_ss + 1, min_fs, max_fs + 1
+ modno, min_ss, max_ss, min_fs, max_fs
) for ((min_fs, max_fs), (min_ss, max_ss)) in delta_method(arr_2d)]
if arr.ndim == 2:
@@ -217,7 +217,7 @@ def to_mask_array(self):
def make_crystfel_bad_regions(self, panels_dict):
modno_to_panels = {}
- for pname, pinfo in panels_dict:
+ for pname, pinfo in panels_dict.items():
modno = crystfel_fmt.panel_modno(pinfo, pname)
modno_to_panels.setdefault(modno, []).append((pname, pinfo))
@@ -303,7 +303,7 @@ def write_crystfel_geom(self, filename, write_mode):
new_bad_regions = self.make_crystfel_bad_regions(geom_dict['panels'])
for i, bad_d in enumerate(new_bad_regions, start=n_area_start):
text_mask.extend([
- f'bad_area{i}/{k} = {v}' for (k, v) in bad_d.items()
+ f'bad_area{i}/{k} = {v}\n' for (k, v) in bad_d.items()
] + ['\n'])
text_write = "".join(text_before + text_mask + text_after)
From 4840021c267eaebd5eba1016e04c724dce56911c Mon Sep 17 00:00:00 2001
From: Thomas Kluyver <thomas@kluyver.me.uk>
Date: Fri, 25 Jun 2021 17:31:04 +0100
Subject: [PATCH 6/9] Give mask objects more helpful reprs
---
extra_geom/mask.py | 8 ++++++++
1 file changed, 8 insertions(+)
diff --git a/extra_geom/mask.py b/extra_geom/mask.py
index 08c4e7d2..9839cad1 100644
--- a/extra_geom/mask.py
+++ b/extra_geom/mask.py
@@ -88,6 +88,9 @@ def __init__(self, modno, start_ss, stop_ss, start_fs, stop_fs):
def _tuple(self):
return self.modno, self.start_ss, self.stop_ss, self.start_fs, self.stop_fs
+ def __repr__(self):
+ return f'RegionRect{self._tuple()}'
+
def __hash__(self):
return hash(self._tuple())
@@ -129,6 +132,11 @@ def __init__(self, shape, regions=()):
self.shape = shape
self.regions = regions # [RegionRect]
+ def __repr__(self):
+ npx = self.to_mask_array().sum()
+ return (f"<MaskRegions for {self.shape} array: "
+ f"{len(self.regions)} covering {npx} pixels>")
+
@classmethod
def from_crystfel_geom(cls, filename):
"""Read mask from a CrystFEL format geometry file.
From 4f8dee87857781556d0486062419e3f9ecce40de Mon Sep 17 00:00:00 2001
From: Thomas Kluyver <thomas@kluyver.me.uk>
Date: Fri, 25 Jun 2021 17:33:59 +0100
Subject: [PATCH 7/9] Clearer repr for MaskRegions
---
extra_geom/mask.py | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/extra_geom/mask.py b/extra_geom/mask.py
index 9839cad1..57d9c070 100644
--- a/extra_geom/mask.py
+++ b/extra_geom/mask.py
@@ -135,7 +135,7 @@ def __init__(self, shape, regions=()):
def __repr__(self):
npx = self.to_mask_array().sum()
return (f"<MaskRegions for {self.shape} array: "
- f"{len(self.regions)} covering {npx} pixels>")
+ f"{len(self.regions)} regions covering {npx} pixels>")
@classmethod
def from_crystfel_geom(cls, filename):
From 7101fc1ad03e5ea06578b18e9a5e967725c36c33 Mon Sep 17 00:00:00 2001
From: Thomas Kluyver <thomas@kluyver.me.uk>
Date: Fri, 25 Jun 2021 18:06:13 +0100
Subject: [PATCH 8/9] Add test for mask -> array and vice versa
---
extra_geom/tests/test_mask.py | 24 ++++++++++++++++++++++++
1 file changed, 24 insertions(+)
diff --git a/extra_geom/tests/test_mask.py b/extra_geom/tests/test_mask.py
index 20317541..d26af5ec 100644
--- a/extra_geom/tests/test_mask.py
+++ b/extra_geom/tests/test_mask.py
@@ -20,3 +20,27 @@ def test_delta_method():
D[0:1, 0:1] = D[2:3, 2:3] = True
exp_res_D = [((0, 1), (0, 1)), ((2, 3), (2, 3))]
assert mask.delta_method(D) == exp_res_D, "Test on separate pixels."
+
+
+def simple_mask():
+ return mask.MaskRegions((2, 4, 5), [
+ mask.RegionRect(None, 0, 1, 0, 2),
+ mask.RegionRect(0, 2, 4, 0, 1),
+ mask.RegionRect(1, 3, 4, 0, 5),
+ ])
+
+def test_roundtrip_array():
+ m1 = simple_mask()
+ arr = m1.to_mask_array()
+ np.testing.assert_array_equal(arr, np.array([
+ [[1, 1, 0, 0, 0],
+ [0, 0, 0, 0, 0],
+ [1, 0, 0, 0, 0],
+ [1, 0, 0, 0, 0]],
+ [[1, 1, 0, 0, 0],
+ [0, 0, 0, 0, 0],
+ [0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1]],
+ ], dtype=np.bool_))
+ m2 = mask.MaskRegions.from_mask_array(arr)
+ assert m2 == m1
From c4b9a3c8e5cbe3d4ec9afe5cdb9b14e5f444c3c5 Mon Sep 17 00:00:00 2001
From: Thomas Kluyver <thomas@kluyver.me.uk>
Date: Thu, 15 Jul 2021 15:30:44 +0100
Subject: [PATCH 9/9] Allow writing geometry object plus mask to .geom file
---
extra_geom/crystfel_fmt.py | 19 +++++++++++-
extra_geom/mask.py | 59 ++++++++++++++++++++++++++------------
2 files changed, 58 insertions(+), 20 deletions(-)
diff --git a/extra_geom/crystfel_fmt.py b/extra_geom/crystfel_fmt.py
index cbfb8af4..dc9bd34f 100644
--- a/extra_geom/crystfel_fmt.py
+++ b/extra_geom/crystfel_fmt.py
@@ -4,6 +4,8 @@
import numpy as np
+from .mask import RegionRect
+
HEADER_TEMPLATE = """\
; {detector} geometry file written by EXtra-geom {version}
; You may need to edit this file to add:
@@ -79,7 +81,8 @@ def frag_to_crystfel(fragment, p, a, ss_slice, fs_slice, dims, pixel_size):
def write_crystfel_geom(self, filename, *,
data_path='/entry_1/instrument_1/detector_1/data',
- mask_path=None, dims=('frame', 'modno', 'ss', 'fs'),
+ mask_path=None, mask_regions=None,
+ dims=('frame', 'modno', 'ss', 'fs'),
nquads=4, adu_per_ev=None, clen=None,
photon_energy=None):
"""Write this geometry to a CrystFEL format (.geom) geometry file.
@@ -116,9 +119,13 @@ def write_crystfel_geom(self, filename, *,
raise ValueError('No frame dimension given')
panel_chunks = []
+ panel_rects = {}
for p, module in enumerate(self.modules):
for a, fragment in enumerate(module):
ss_slice, fs_slice = self._tile_slice(a)
+ panel_rects[f'p{p}a{a}'] = RegionRect(
+ p, ss_slice.start, ss_slice.stop, fs_slice.start, fs_slice.stop
+ )
if 'modno' not in dims:
# If we don't have a modno dimension, assume modules are
# concatenated along the slow-scan dim, e.g. AGIPD (8192, 128)
@@ -139,6 +146,15 @@ def write_crystfel_geom(self, filename, *,
paths['mask'] = mask_path
path_str = '\n'.join('{} = {} ;'.format(i, j) for i, j in paths.items())
+ mask_lines = []
+ if mask_regions is not None:
+ for i, bad_d in enumerate(mask_regions.make_crystfel_bad_regions(
+ panel_rects, modules_stacked=('modno' in dims)
+ )):
+ mask_lines.extend([
+ f'bad_area{i}/{k} = {v}' for (k, v) in bad_d.items()
+ ] + [''])
+
with open(filename, 'w') as f:
f.write(HEADER_TEMPLATE.format(
detector=self.detector_type_name,
@@ -152,6 +168,7 @@ def write_crystfel_geom(self, filename, *,
))
rigid_groups = get_rigid_groups(self, nquads=nquads)
f.write(rigid_groups)
+ f.write('\n'.join(mask_lines))
for chunk in panel_chunks:
f.write(chunk)
diff --git a/extra_geom/mask.py b/extra_geom/mask.py
index 57d9c070..ac9c8f0e 100644
--- a/extra_geom/mask.py
+++ b/extra_geom/mask.py
@@ -97,6 +97,13 @@ def __hash__(self):
def __eq__(self, other):
return isinstance(other, RegionRect) and self._tuple() == other._tuple()
+ def replace(self, **kwargs):
+ new_kwargs = {
+ k: kwargs.get(k, getattr(self, k))
+ for k in ('modno', 'start_ss', 'stop_ss', 'start_fs', 'stop_fs')
+ }
+ return RegionRect(**new_kwargs)
+
@property
def array_slice(self):
"""Get a tuple to use for slicing an array"""
@@ -223,13 +230,18 @@ def to_mask_array(self):
return res_mask
- def make_crystfel_bad_regions(self, panels_dict):
- modno_to_panels = {}
- for pname, pinfo in panels_dict.items():
- modno = crystfel_fmt.panel_modno(pinfo, pname)
- modno_to_panels.setdefault(modno, []).append((pname, pinfo))
+ def make_crystfel_bad_regions(self, panel_rects, modules_stacked=True):
+ panel_rects_by_modno = {i: {} for i in range(self.shape[0])}
+ for pname, region in panel_rects.items():
+ panel_rects_by_modno[region.modno][pname] = region
def to_dict(region: RegionRect):
+ if (region.modno is not None) and not modules_stacked:
+ mod_offset = region.modno * self.shape[1]
+ region = region.replace(
+ start_ss=region.start_ss + mod_offset,
+ stop_ss=region.stop_ss + mod_offset,
+ )
return {
'min_ss': region.start_ss, 'max_ss': region.stop_ss - 1,
'min_fs': region.start_fs, 'max_fs': region.stop_fs - 1,
@@ -237,21 +249,25 @@ def to_dict(region: RegionRect):
res = []
for mask_region in self.regions:
- if mask_region.modno is None:
- res.append(to_dict(mask_region)) # Mask for all panels
+ if modules_stacked:
+ if mask_region.modno is None:
+ yield to_dict(mask_region) # Mask for all panels
+ else:
+ module_panels = panel_rects_by_modno[mask_region.modno]
+ for pname, panel_region in module_panels.items():
+ overlap, matching = mask_region.intersection(panel_region)
+ if overlap:
+ d = to_dict(matching)
+ d['panel'] = pname
+ res.append(d)
+
else:
- for pname, pinfo in modno_to_panels[mask_region.modno]:
- panel_region = RegionRect(
- mask_region.modno,
- pinfo['min_ss'], pinfo['max_ss'] + 1,
- pinfo['min_fs'], pinfo['max_fs'] + 1,
- )
- overlap, matching = mask_region.intersection(panel_region)
- if overlap:
- d = to_dict(matching)
- d['panel'] = pname
- res.append(d)
+ if mask_region.modno is None:
+ for i in range(self.shape[0]):
+ yield to_dict(mask_region.replace(modno=i))
+ else:
+ yield to_dict(mask_region)
return res
@@ -308,7 +324,12 @@ def write_crystfel_geom(self, filename, write_mode):
geom_dict = load_crystfel_geometry(filename)
- new_bad_regions = self.make_crystfel_bad_regions(geom_dict['panels'])
+ panel_rects = {pname: RegionRect(
+ crystfel_fmt.panel_modno(pinfo, pname),
+ pinfo['min_ss'], pinfo['max_ss'] + 1,
+ pinfo['min_fs'], pinfo['max_fs'] + 1,
+ ) for (pname, pinfo) in geom_dict['panels'].items()}
+ new_bad_regions = self.make_crystfel_bad_regions(panel_rects)
for i, bad_d in enumerate(new_bad_regions, start=n_area_start):
text_mask.extend([
f'bad_area{i}/{k} = {v}\n' for (k, v) in bad_d.items()