Merge pull request #2932 from cta-observatory/camera-pixel-grid-type

Implement square pixels on a hexgrid

Author: maxnoe

docs/changes/2932.feature.rst

@@ -0,0 +1,2 @@
+Add a new attribute ``grid_type`` to the ``CameraGeometry`` to
+add the possibility of handling cameras with square pixels using hexagonal grids.

pyproject.toml

@@ -176,6 +176,7 @@ filterwarnings = [
     "ignore:'setParseAction':DeprecationWarning:matplotlib",
     "ignore:'endQuoteChar':DeprecationWarning:matplotlib",
     "ignore:'unquoteResults':DeprecationWarning:matplotlib",
+    "ignore:'enablePackrat':DeprecationWarning:matplotlib",
     "ignore:'parseAll':DeprecationWarning:pyparsing",
 ]
 norecursedirs = [

src/ctapipe/instrument/__init__.py

@@ -1,5 +1,11 @@
 from .atmosphere import get_atmosphere_profile_functions
-from .camera import CameraDescription, CameraGeometry, CameraReadout, PixelShape
+from .camera import (
+    CameraDescription,
+    CameraGeometry,
+    CameraReadout,
+    PixelGridType,
+    PixelShape,
+)
 from .guess import guess_telescope
 from .optics import (
     ComaPSFModel,
@@ -21,6 +27,7 @@
     "get_atmosphere_profile_functions",
     "guess_telescope",
     "OpticsDescription",
+    "PixelGridType",
     "PixelShape",
     "SubarrayDescription",
     "TelescopeDescription",

src/ctapipe/instrument/camera/__init__.py

@@ -1,5 +1,10 @@
 from .description import CameraDescription  # noqa: F401
-from .geometry import CameraGeometry, PixelShape, UnknownPixelShapeWarning  # noqa: F401
+from .geometry import (
+    CameraGeometry,  # noqa: F401
+    PixelGridType,  # noqa: F401
+    PixelShape,  # noqa: F401
+    UnknownPixelShapeWarning,  # noqa: F401
+)
 from .readout import CameraReadout  # noqa: F401
 
 # commented out due to sphinx issue with classes being defined in 3 places

src/ctapipe/instrument/camera/geometry.py

@@ -14,7 +14,7 @@
 from astropy.table import Table
 from astropy.utils import lazyproperty
 from scipy.sparse import csr_array, lil_array
-from scipy.spatial import cKDTree
+from scipy.spatial import KDTree
 
 from ctapipe.coordinates import CameraFrame, get_representation_component_names
 from ctapipe.utils import get_table_dataset
@@ -36,6 +36,18 @@ def _distance(x1, y1, x2, y2):
     return np.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
 
 
+@unique
+class PixelGridType(Enum):
+    """Grid type for pixel coordinate grids."""
+
+    #: Regular square grid with minimal gaps
+    REGULAR_SQUARE = "square"
+    #: Regular hexagonal grid
+    REGULAR_HEX = "hex"
+    #: Everything else
+    IRREGULAR = "irregular"
+
+
 @unique
 class PixelShape(Enum):
     """Supported Pixel Shapes Enum"""
@@ -100,6 +112,10 @@ class CameraGeometry:
         position of each pixel (y-coordinate)
     pix_area : u.Quantity
         surface area of each pixel
+    grid_type : None | PixelGridType
+        The grid type of the pixel grid. By default, square pixels
+        are assumed to be on a regular square grid and circular and hexagonal
+        pixels are assumed to be on a regular hex grid.
     pix_type : PixelShape
         either 'rectangular' or 'hexagonal'
     pix_rotation : u.Quantity[angle]
@@ -118,8 +134,8 @@ class CameraGeometry:
         Frame in which the pixel coordinates are defined (after applying cam_rotation)
     """
 
-    CURRENT_TAB_VERSION = "2.0"
-    SUPPORTED_TAB_VERSIONS = {"1.0", "1", "1.1", "2.0"}
+    CURRENT_TAB_VERSION = "2.1"
+    SUPPORTED_TAB_VERSIONS = {"1.0", "1", "1.1", "2.0", "2.1"}
 
     def __init__(
         self,
@@ -128,7 +144,8 @@ def __init__(
         pix_x,
         pix_y,
         pix_area,
-        pix_type,
+        pix_type: PixelShape,
+        grid_type: None | PixelGridType | str = None,
         pix_rotation=0 * u.deg,
         cam_rotation=0 * u.deg,
         neighbors=None,
@@ -140,52 +157,19 @@ def __init__(
         self.n_pixels = len(pix_id)
         self.unit = pix_x.unit
         self.pix_id = np.array(pix_id)
-
-        if _validate:
-            self.pix_id = np.array(self.pix_id)
-
-            if self.pix_id.ndim != 1:
-                raise ValueError(
-                    f"Pixel coordinates must be 1 dimensional, got {pix_id.ndim}"
-                )
-
-            shape = (self.n_pixels,)
-
-            if pix_x.shape != shape:
-                raise ValueError(
-                    f"pix_x has wrong shape: {pix_x.shape}, expected {shape}"
-                )
-            if pix_y.shape != shape:
-                raise ValueError(
-                    f"pix_y has wrong shape: {pix_y.shape}, expected {shape}"
-                )
-            if pix_area.shape != shape:
-                raise ValueError(
-                    f"pix_area has wrong shape: {pix_area.shape}, expected {shape}"
-                )
-
-            if isinstance(pix_type, str):
-                pix_type = PixelShape.from_string(pix_type)
-            elif not isinstance(pix_type, PixelShape):
-                raise TypeError(
-                    f"pix_type must be a PixelShape or the name of a PixelShape, got {pix_type}"
-                )
-
-            if not isinstance(pix_rotation, Angle):
-                pix_rotation = Angle(pix_rotation)
-
-            if not isinstance(cam_rotation, Angle):
-                cam_rotation = Angle(cam_rotation)
-
         self.pix_x = pix_x
         self.pix_y = pix_y.to(self.unit)
         self.pix_area = pix_area.to(self.unit**2)
         self.pix_type = pix_type
+        self.grid_type = grid_type
         self.pix_rotation = pix_rotation
         self.cam_rotation = cam_rotation
         self._neighbors = neighbors
         self.frame = frame
 
+        if _validate:
+            self._validate()
+
         if neighbors is not None:
             if isinstance(neighbors, list):
                 lil = lil_array((self.n_pixels, self.n_pixels), dtype=bool)
@@ -201,6 +185,59 @@ def __init__(
             # cache border pixel mask per instance
         self._border_cache = {}
 
+    def _validate(self):
+        self.pix_id = np.array(self.pix_id)
+
+        if self.pix_id.ndim != 1:
+            raise ValueError(
+                f"Pixel coordinates must be 1 dimensional, got {self.pix_id.ndim}"
+            )
+
+        shape = (self.n_pixels,)
+
+        if self.pix_x.shape != shape:
+            raise ValueError(
+                f"pix_x has wrong shape: {self.pix_x.shape}, expected {shape}"
+            )
+        if self.pix_y.shape != shape:
+            raise ValueError(
+                f"pix_y has wrong shape: {self.pix_y.shape}, expected {shape}"
+            )
+        if self.pix_area.shape != shape:
+            raise ValueError(
+                f"pix_area has wrong shape: {self.pix_area.shape}, expected {shape}"
+            )
+
+        if isinstance(self.pix_type, str):
+            self.pix_type = PixelShape.from_string(self.pix_type)
+        elif not isinstance(self.pix_type, PixelShape):
+            raise TypeError(
+                f"pix_type must be a PixelShape or the name of a PixelShape, got {self.pix_type}"
+            )
+
+        # if grid_type is not given, deduce grid type from pix type assuming regular grids
+        self.grid_type = self._get_grid_type(self.grid_type, self.pix_type)
+
+        if not isinstance(self.pix_rotation, Angle):
+            self.pix_rotation = Angle(self.pix_rotation)
+
+        if not isinstance(self.cam_rotation, Angle):
+            self.cam_rotation = Angle(self.cam_rotation)
+
+    @staticmethod
+    def _get_grid_type(grid_type: PixelGridType | str | None, pixel_shape: PixelShape):
+        if isinstance(grid_type, PixelGridType):
+            return grid_type
+
+        if grid_type is None:
+            # backwards compatibility: assume square grid for square pixels, hexgrid for anything else
+            if pixel_shape is PixelShape.SQUARE:
+                return PixelGridType.REGULAR_SQUARE
+
+            return PixelGridType.REGULAR_HEX
+
+        return PixelGridType(grid_type)
+
     def __eq__(self, other):
         if not isinstance(other, CameraGeometry):
             return NotImplemented
@@ -323,6 +360,7 @@ def transform_to(self, frame: BaseCoordinateFrame):
             pix_y=trans_y,
             pix_area=pix_area,
             pix_type=cam.pix_type,
+            grid_type=cam.grid_type,
             pix_rotation=pix_rotation,
             cam_rotation=cam_rotation,
             neighbors=cam._neighbors,
@@ -352,6 +390,7 @@ def __getitem__(self, slice_):
             pix_y=self.pix_y[slice_],
             pix_area=self.pix_area[slice_],
             pix_type=self.pix_type,
+            grid_type=self.grid_type,
             pix_rotation=self.pix_rotation,
             cam_rotation=self.cam_rotation,
             neighbors=None,
@@ -425,7 +464,7 @@ def _kdtree(self):
         """
 
         pixel_centers = np.column_stack([self.pix_x.value, self.pix_y.value])
-        return cKDTree(pixel_centers)
+        return KDTree(pixel_centers)
 
     @lazyproperty
     def _all_pixel_areas_equal(self):
@@ -621,6 +660,7 @@ def to_table(self):
             names=["pix_id", "pix_x", "pix_y", "pix_area"],
             meta=dict(
                 PIX_TYPE=self.pix_type.value,
+                GRID=self.grid_type.value,
                 TAB_TYPE="ctapipe.instrument.CameraGeometry",
                 TAB_VER=self.CURRENT_TAB_VERSION,
                 CAM_ID=self.name,
@@ -669,6 +709,7 @@ def from_table(cls, url_or_table, **kwargs):
             pix_y=tab["pix_y"].quantity,
             pix_area=tab["pix_area"].quantity,
             pix_type=tab.meta["PIX_TYPE"],
+            grid_type=tab.meta.get("GRID"),
             pix_rotation=Angle(tab.meta["PIX_ROT"], u.deg),
             cam_rotation=Angle(tab.meta["CAM_ROT"], u.deg),
         )
@@ -722,15 +763,17 @@ def calc_pixel_neighbors(self, diagonal=False):
         if self.n_pixels <= 1:
             return csr_array(np.ones((self.n_pixels, self.n_pixels), dtype=bool))
 
-        # assume circle pixels are also on a hex grid
-        if self.pix_type in (PixelShape.HEXAGON, PixelShape.CIRCLE):
+        if self.grid_type is PixelGridType.REGULAR_HEX:
             max_neighbors = 6
             # on a hexgrid, the closest pixel in the second circle is
             # the diameter of the hexagon plus the inradius away
             # in units of the diameter, this is 1 + np.sqrt(3) / 4 = 1.433
             radius = 1.4
             norm = 2  # use L2 norm for hex
-        else:
+            # for square pixels on a hexgrid, the grid is stretched
+            if self.pix_type is PixelShape.SQUARE:
+                radius *= np.sqrt(5) / 2
+        elif self.grid_type is PixelGridType.REGULAR_SQUARE:
             # if diagonal should count as neighbor, we
             # need to find at most 8 neighbors with a max L2 distance
             # < than 2 * the pixel size, else 4 neighbors with max L1 distance
@@ -745,6 +788,10 @@ def calc_pixel_neighbors(self, diagonal=False):
                 max_neighbors = 4
                 radius = 1.5
                 norm = 1
+        else:
+            raise ValueError(
+                "Automatic computation of pixel neighbors only implemented for regular square and hex grids"
+            )
 
         distances, neighbor_candidates = self._kdtree.query(
             self._kdtree.data, k=max_neighbors + 1, p=norm

src/ctapipe/instrument/camera/tests/test_geometry.py

@@ -407,3 +407,19 @@ def test_empty(prod5_lst):
     empty = geometry[mask]
 
     assert empty.neighbor_matrix.shape == (0, 0)
+
+
+def test_squre_pixels_on_hexgrid(geometry_hexgrid_square_pixels):
+    geom = geometry_hexgrid_square_pixels
+
+    border = geom.get_border_pixel_mask()
+    # 16 by 16 pixels means 60 border pixels 4 * 14 edge pixels + 4 corrner pixels
+    assert np.count_nonzero(border) == 4 * 14 + 4
+
+    n_neighbors = np.count_nonzero(geom.neighbor_matrix, axis=0)
+    np.testing.assert_array_equal(n_neighbors[~border], 6)
+
+    table = geom.to_table()
+    assert table.meta["GRID"] == "hex"
+    geom_from_table = CameraGeometry.from_table(table)
+    assert geom_from_table == geom

src/ctapipe/instrument/conftest.py

@@ -5,6 +5,8 @@
 import os
 import shutil
 
+import astropy.units as u
+import numpy as np
 import pytest
 
 from ctapipe.utils.filelock import FileLock
@@ -42,3 +44,37 @@ def svc_path(instrument_dir):
         del os.environ["CTAPIPE_SVC_PATH"]
     else:
         os.environ["CTAPIPE_SVC_PATH"] = before
+
+
+@pytest.fixture()
+def geometry_hexgrid_square_pixels():
+    """A camera with square pixels on a hexagonal grid"""
+    from ctapipe.coordinates import CameraFrame
+    from ctapipe.instrument import CameraGeometry, PixelGridType, PixelShape
+
+    size = 22
+    pix_id = np.arange(256)
+
+    coords = np.arange(-7.5 * size, 7.6 * size, size)
+    pix_x, pix_y = np.meshgrid(coords, coords)
+
+    # offset every second row by half the pixel size
+    pix_x[::2] += 0.5 * size
+    pix_x = pix_x.ravel() * u.mm
+    pix_y = pix_y.ravel() * u.mm
+
+    # introduce some gaps
+    pix_area = (size / 1.05) ** 2
+    pix_area = np.full(len(pix_id), pix_area) * u.mm**2
+
+    geom = CameraGeometry(
+        pix_id=pix_id,
+        pix_x=pix_x,
+        pix_y=pix_y,
+        pix_area=pix_area,
+        pix_type=PixelShape.SQUARE,
+        grid_type=PixelGridType.REGULAR_HEX,
+        name="HEXSQUARECAM",
+        frame=CameraFrame(focal_length=10 * u.m),
+    )
+    return geom