Merge pull request #111 from fact-project/neighbors

Add function to return pixel neighborhood matrix

fact/VERSION

@@ -1 +1 @@
-0.20.1
+0.21.0

fact/instrument/camera.py

@@ -2,10 +2,12 @@
 import numpy as np
 from functools import lru_cache
 import pandas as pd
+from scipy.sparse import csr_matrix
+from scipy.spatial import cKDTree
 
 from .constants import (
     FOCAL_LENGTH_MM, PINCUSHION_DISTORTION_SLOPE,
-    PIXEL_SPACING_MM, FOV_PER_PIXEL_DEG
+    PIXEL_SPACING_MM, FOV_PER_PIXEL_DEG, N_PIXEL
 )
 
 
@@ -172,8 +174,10 @@ def combine_bias_patch_current_to_trigger_patch_current(bias_patch_currents):
     fivers = pi[pi.bias_patch_size == 5].sort_values('trigger_patch_id')
 
     b_c = bias_patch_currents  # just to shorten the name
-    t_c = b_c[fourers.bias_patch_id.values] * 4/9 + b_c[fivers.bias_patch_id.values] * 5/9
-
+    t_c = (
+        b_c[fourers.bias_patch_id.values] * 4 / 9
+        + b_c[fivers.bias_patch_id.values] * 5 / 9
+    )
     trigger_patch_currents = t_c  # unshorten the name
     return trigger_patch_currents
 
@@ -190,3 +194,39 @@ def take_apart_trigger_values_for_bias_patches(trigger_rates):
     pi = patch_indices().sort_values('bias_patch_id')
 
     return trigger_rates[pi.trigger_patch_id.values]
+
+
+@lru_cache(maxsize=1)
+def get_neighbor_matrix():
+    '''
+    Returns a sparse boolean neighbor matrix with n[chid, other_chid] = is neighbor.
+    '''
+    xy = get_pixel_dataframe().loc[:, ['x', 'y']].values
+    tree = cKDTree(xy)
+    neighbors = tree.query_ball_tree(tree, r=10)
+
+    n_neighbors = [len(n) for n in neighbors]
+    col = np.repeat(np.arange(N_PIXEL), n_neighbors)
+    row = [pix for n in neighbors for pix in n]
+    data = np.ones(len(row))
+    m = csr_matrix((data, (row, col)), shape=(N_PIXEL, N_PIXEL), dtype=bool)
+    m.setdiag(False)
+    return m
+
+
+@lru_cache(maxsize=1)
+def get_num_neighbors():
+    '''
+    Return a numpy array with the number of neighbors for each pixel in chid order
+    '''
+    return get_neighbor_matrix().sum(axis=0).A1
+
+
+@lru_cache()
+def get_border_pixel_mask(width=1):
+    if width == 1:
+        return get_num_neighbors() < 6
+
+    n = get_neighbor_matrix().todense().A
+
+    return (n & get_border_pixel_mask(width - 1)).any(axis=1)

fact/instrument/constants.py

@@ -5,6 +5,7 @@
 from astropy.coordinates import EarthLocation
 import astropy.units as u
 
+N_PIXEL = 1440
 #: The inner diameter of the hexagonal pixels in mm.
 #: This is also the grid constant of the hex grid.
 PIXEL_SPACING_MM = 9.5

tests/test_camera.py

@@ -0,0 +1,28 @@
+import numpy as np
+
+
+def test_neighbors():
+    from fact.instrument.camera import get_neighbor_matrix
+
+    neighbors = get_neighbor_matrix()
+
+    assert neighbors[1144, 259]
+    assert neighbors[1144, 1143]
+    assert neighbors[1144, 1146]
+    assert neighbors[1144, 1147]
+    assert neighbors[1144, 287]
+    assert neighbors[1144, 284]
+    assert not neighbors[1144, 256]
+    assert not neighbors[1144, 281]
+
+    assert np.all(neighbors.diagonal() == 0)
+
+
+def test_n_neighbors():
+    from fact.instrument.camera import get_num_neighbors
+
+    n_neighbors = get_num_neighbors()
+
+    assert n_neighbors[54] == 3
+    assert n_neighbors[86] == 3
+    assert n_neighbors[81] == 4