Fixing the toy image generator (#790)

* Add test for toy image generator

* Fix toymodel generator and test event

* The evaluation of the 2d gaussian was missing the normalisation
for the pixel area
* Length and width were not squared
* Adapted test image to show a more realistic shower image

* Remove print

* Adapt hillas example

* Adapt shower image size to camera geometry

* Use larger tolerance

* Adapt test image to camera geometry

* Use itertools for test loop

* Make hillas test agnostic to 180deg rotations

* Adapt examples to new toy generator

ctapipe/image/tests/test_geometry_converter.py

@@ -12,19 +12,33 @@
 cam_ids = CameraGeometry.get_known_camera_names()
 
 
+def create_mock_image(geom):
+    '''
+    creates a mock image, which parameters are adapted to the camera size
+    '''
+
+    camera_r = np.max(np.sqrt(geom.pix_x**2 + geom.pix_y**2))
+    model = generate_2d_shower_model(
+        centroid=(0.3 * camera_r.value, 0),
+        width=0.03 * camera_r.value,
+        length=0.10 * camera_r.value,
+        psi="25d"
+    )
+
+    _, image, _ = make_toymodel_shower_image(
+        geom, model.pdf,
+        intensity=0.5 * geom.n_pixels,
+        nsb_level_pe=3,
+    )
+    return image
+
+
 @pytest.mark.parametrize("rot", [3, ])
 @pytest.mark.parametrize("cam_id", cam_ids)
 def test_convert_geometry(cam_id, rot):
 
     geom = CameraGeometry.from_name(cam_id)
-
-    model = generate_2d_shower_model(centroid=(0.4, 0), width=0.01, length=0.03,
-                                     psi="25d")
-
-    _, image, _ = make_toymodel_shower_image(geom, model.pdf,
-                                             intensity=50,
-                                             nsb_level_pe=100)
-
+    image = create_mock_image(geom)
     hillas_0 = hillas_parameters(geom, image)
 
     if geom.pix_type == 'hexagonal':
@@ -70,14 +84,7 @@ def test_convert_geometry_mock(cam_id, rot):
     """
 
     geom = CameraGeometry.from_name(cam_id)
-
-    model = generate_2d_shower_model(centroid=(0.4, 0), width=0.01, length=0.03,
-                                     psi="25d")
-
-    _, image, _ = make_toymodel_shower_image(geom, model.pdf,
-                                             intensity=50,
-                                             nsb_level_pe=100)
-
+    image = create_mock_image(geom)
     hillas_0 = hillas_parameters(geom, image)
 
     if geom.pix_type == 'hexagonal':

ctapipe/image/tests/test_hillas.py

@@ -9,6 +9,7 @@
 from numpy.random import seed
 from numpy.ma import masked_array
 import pytest
+from itertools import combinations
 
 methods = (
     hillas_parameters_1,
@@ -28,18 +29,21 @@ def create_sample_image(psi='-30d'):
     geom = CameraGeometry.from_name("LSTCam")
 
     # make a toymodel shower model
-    model = toymodel.generate_2d_shower_model(centroid=(0.2, 0.3),
-                                              width=0.001, length=0.01,
-                                              psi=psi)
+    model = toymodel.generate_2d_shower_model(
+        centroid=(0.2, 0.3),
+        width=0.05, length=0.15,
+        psi=psi,
+    )
 
     # generate toymodel image in camera for this shower model.
-    image, signal, noise = toymodel.make_toymodel_shower_image(geom, model.pdf,
-                                                               intensity=50,
-                                                               nsb_level_pe=100)
+    image, signal, noise = toymodel.make_toymodel_shower_image(
+        geom, model.pdf,
+        intensity=1500,
+        nsb_level_pe=3,
+    )
 
     # denoise the image, so we can calculate hillas params
-    clean_mask = tailcuts_clean(geom, image, 10,
-                                5)  # pedvars = 1 and core and boundary
+    clean_mask = tailcuts_clean(geom, image, 10, 5)
 
     return geom, image, clean_mask
 
@@ -86,18 +90,24 @@ def test_hillas():
             for i, method in enumerate(methods, start=1)
         }
 
+        for result in results.values():
+            if result.psi < -90 * u.deg:
+                result.psi += 180 * u.deg
+                result.skewness *= -1
+            elif result.psi > 90 * u.deg:
+                result.psi -= 180 * u.deg
+                result.skewness *= -1
+
         # compare each method's output
-        for aa in results:
-            for bb in results:
-                if aa is not bb:
-                    print("comparing {} to {}".format(aa, bb))
-                    compare_result(results[aa].length, results[bb].length)
-                    compare_result(results[aa].width, results[bb].width)
-                    compare_result(results[aa].r, results[bb].r)
-                    compare_result(results[aa].phi.deg, results[bb].phi.deg)
-                    compare_result(results[aa].psi.deg, results[bb].psi.deg)
-                    compare_result(results[aa].skewness, results[bb].skewness)
-                    # compare_result(results[aa].kurtosis, results[bb].kurtosis)
+        for aa, bb in combinations(results, 2):
+            print("comparing {} to {}".format(aa, bb))
+            compare_result(results[aa].length, results[bb].length)
+            compare_result(results[aa].width, results[bb].width)
+            compare_result(results[aa].r, results[bb].r)
+            compare_result(results[aa].phi.deg, results[bb].phi.deg)
+            compare_result(results[aa].psi.deg, results[bb].psi.deg)
+            compare_result(results[aa].skewness, results[bb].skewness)
+            # compare_result(results[aa].kurtosis, results[bb].kurtosis)
 
 
 def test_hillas_masked():

ctapipe/image/tests/test_toy.py

@@ -0,0 +1,41 @@
+import numpy as np
+from ctapipe.instrument import CameraGeometry
+from pytest import approx
+from scipy.stats import poisson
+
+
+def test_intensity():
+    from .. import toymodel
+
+    np.random.seed(0)
+
+    geom = CameraGeometry.from_name('LSTCam')
+
+    width = 0.05
+    length = 0.15
+    intensity = 50
+
+    # make a toymodel shower model
+    model = toymodel.generate_2d_shower_model(
+        centroid=(0.2, 0.3),
+        width=width, length=length,
+        psi='30d',
+    )
+
+    image, signal, noise = toymodel.make_toymodel_shower_image(
+        geom, model.pdf, intensity=intensity, nsb_level_pe=5,
+    )
+
+    # test if signal reproduces given cog values
+    assert np.average(geom.pix_x.value, weights=signal) == approx(0.2, rel=0.15)
+    assert np.average(geom.pix_y.value, weights=signal) == approx(0.3, rel=0.15)
+
+    # test if signal reproduces given width/length values
+    cov = np.cov(geom.pix_x.value, geom.pix_y.value, aweights=signal)
+    eigvals, eigvecs = np.linalg.eigh(cov)
+
+    assert np.sqrt(eigvals[0]) == approx(width, rel=0.15)
+    assert np.sqrt(eigvals[1]) == approx(length, rel=0.15)
+
+    # test if total intensity is inside in 99 percent confidence interval
+    assert poisson(intensity).ppf(0.05) <= signal.sum() <= poisson(intensity).ppf(0.95)

ctapipe/image/toymodel.py

@@ -1,6 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """
-Utilities to generate toymodel (fake) reconstruction inputs for testing 
+Utilities to generate toymodel (fake) reconstruction inputs for testing
 purposes.
 
 Example:
@@ -9,12 +9,12 @@
 
     >>> from instrument import CameraGeometry
     >>> geom = CameraGeometry.make_rectangular(20,20)
-    >>> showermodel = generate_2d_shower_model(centroid=[0.25, 0.0], 
+    >>> showermodel = generate_2d_shower_model(centroid=[0.25, 0.0],
     length=0.1,width=0.02, psi='40d')
     >>> image, signal, noise = make_toymodel_shower_image(geom, showermodel.pdf)
     >>> print(image.shape)
     (400,)
-                                             
+
 .. plot:: image/image_example.py
     :include-source:
 
@@ -53,7 +53,7 @@ def generate_2d_shower_model(centroid, width, length, psi):
     a `scipy.stats` object
 
     """
-    aligned_covariance = np.array([[length, 0], [0, width]])
+    aligned_covariance = np.array([[length**2, 0], [0, width**2]])
     # rotate by psi angle: C' = R C R+
     rotation = linalg.rotation_matrix_2d(psi)
     rotated_covariance = rotation.dot(aligned_covariance).dot(rotation.T)
@@ -69,9 +69,9 @@ def make_toymodel_shower_image(geom, showerpdf, intensity=50, nsb_level_pe=50):
     Parameters
     ----------
     geom : `ctapipe.instrument.CameraGeometry`
-        camera geometry object 
+        camera geometry object
     showerpdf : func
-        PDF function for the shower to generate in the camera, e.g. from a 
+        PDF function for the shower to generate in the camera, e.g. from a
     intensity : int
         factor to multiply the model by to get photo-electrons
     nsb_level_pe : type
@@ -84,10 +84,11 @@ def make_toymodel_shower_image(geom, showerpdf, intensity=50, nsb_level_pe=50):
 
     """
     pos = np.empty(geom.pix_x.shape + (2,))
-    pos[..., 0] = geom.pix_x.value
-    pos[..., 1] = geom.pix_y.value
+    pos[:, 0] = geom.pix_x.value
+    pos[:, 1] = geom.pix_y.value
+
+    model_counts = showerpdf(pos) * intensity * geom.pix_area.value
 
-    model_counts = (showerpdf(pos) * intensity).astype(np.int32)
     signal = np.random.poisson(model_counts)
     noise = np.random.poisson(nsb_level_pe, size=signal.shape)
     image = (signal + noise) - np.mean(noise)

examples/camera_animation.py

@@ -22,26 +22,23 @@
 
     # load the camera
     tel = TelescopeDescription.from_name("SST-1M", "DigiCam")
-    print(tel, tel.optics.equivalent_focal_length)
     geom = tel.camera
 
-    # poor-man's coordinate transform from telscope to camera frame (it's
-    # better to use ctapipe.coordiantes when they are stable)
-    scale = tel.optics.equivalent_focal_length.to(geom.pix_x.unit).value
-    fov = np.deg2rad(4.0)
-    maxwid = np.deg2rad(0.01)
-    maxlen = np.deg2rad(0.03)
+    fov = 0.3
+    maxwid = 0.05
+    maxlen = 0.1
 
     disp = CameraDisplay(geom, ax=ax)
-    disp.cmap = plt.cm.terrain
+    disp.cmap = 'inferno'
     disp.add_colorbar(ax=ax)
 
     def update(frame):
-        centroid = np.random.uniform(-fov, fov, size=2) * scale
-        width = np.random.uniform(0, maxwid) * scale
-        length = np.random.uniform(0, maxlen) * scale + width
-        angle = np.random.uniform(0, 360)
-        intens = np.random.exponential(2) * 50
+        centroid = np.random.uniform(-fov, fov, size=2)
+        width = np.random.uniform(0.01, maxwid)
+        length = np.random.uniform(width, maxlen)
+        angle = np.random.uniform(0, 180)
+        intens = width * length * (5e4 + 1e5 * np.random.exponential(2))
+
         model = toymodel.generate_2d_shower_model(
             centroid=centroid,
             width=width,
@@ -52,10 +49,9 @@ def update(frame):
             geom,
             model.pdf,
             intensity=intens,
-            nsb_level_pe=5000,
+            nsb_level_pe=5,
         )
-        image /= image.max()
         disp.image = image
 
-    anim = FuncAnimation(fig, update, interval=250)
+    anim = FuncAnimation(fig, update, interval=500)
     plt.show()

examples/camera_display.py

@@ -10,36 +10,25 @@
 from ctapipe.visualization import CameraDisplay
 
 
-def draw_neighbors(geom, pixel_index, color='r', **kwargs):
-    """Draw lines between a pixel and its neighbors"""
-
-    neigh = geom.neighbors[pixel_index]  # neighbor indices (not pixel ids)
-    x, y = geom.pix_x[pixel_index].value, geom.pix_y[pixel_index].value
-    for nn in neigh:
-        nx, ny = geom.pix_x[nn].value, geom.pix_y[nn].value
-        plt.plot([x, nx], [y, ny], color=color, **kwargs)
-
-
 if __name__ == '__main__':
 
     # Load the camera
     geom = CameraGeometry.from_name("LSTCam")
     disp = CameraDisplay(geom)
-    disp.set_limits_minmax(0, 300)
     disp.add_colorbar()
 
     # Create a fake camera image to display:
     model = toymodel.generate_2d_shower_model(
-        centroid=(0.2, 0.0), width=0.01, length=0.1, psi='35d'
+        centroid=(0.2, 0.0), width=0.05, length=0.15, psi='35d'
     )
 
     image, sig, bg = toymodel.make_toymodel_shower_image(
-        geom, model.pdf, intensity=50, nsb_level_pe=1000
+        geom, model.pdf, intensity=1500, nsb_level_pe=2
     )
 
     # Apply image cleaning
     cleanmask = tailcuts_clean(
-        geom, image, picture_thresh=200, boundary_thresh=100
+        geom, image, picture_thresh=10, boundary_thresh=5
     )
     clean = image.copy()
     clean[~cleanmask] = 0.0
@@ -50,15 +39,8 @@ def draw_neighbors(geom, pixel_index, color='r', **kwargs):
 
     # Show the camera image and overlay Hillas ellipse and clean pixels
     disp.image = image
-    disp.cmap = 'PuOr'
-    disp.highlight_pixels(cleanmask, color='black')
-    disp.overlay_moments(hillas, color='cyan', linewidth=3)
-
-    # Draw the neighbors of pixel 100 in red, and the neighbor-neighbors in
-    # green
-    for ii in geom.neighbors[130]:
-        draw_neighbors(geom, ii, color='green')
-
-    draw_neighbors(geom, 130, color='cyan', lw=2)
+    disp.cmap = 'inferno'
+    disp.highlight_pixels(cleanmask, color='crimson')
+    disp.overlay_moments(hillas, color='cyan', linewidth=1)
 
     plt.show()

examples/camera_display_multi.py

@@ -19,7 +19,7 @@ def draw_several_cams(geom, ncams=4):
 
     cmaps = ['jet', 'afmhot', 'terrain', 'autumn']
     fig, axs = plt.subplots(
-        1, ncams, figsize=(15, 4), sharey=True, sharex=True
+        1, ncams, figsize=(15, 4),
     )
 
     for ii in range(ncams):
@@ -32,16 +32,16 @@ def draw_several_cams(geom, ncams=4):
 
         model = toymodel.generate_2d_shower_model(
             centroid=(0.2 - ii * 0.1, -ii * 0.05),
-            width=0.005 + 0.001 * ii,
-            length=0.1 + 0.05 * ii,
+            width=0.05 + 0.001 * ii,
+            length=0.15 + 0.05 * ii,
             psi=ii * 20 * u.deg,
         )
 
         image, sig, bg = toymodel.make_toymodel_shower_image(
             geom,
             model.pdf,
-            intensity=50,
-            nsb_level_pe=1000,
+            intensity=1500,
+            nsb_level_pe=5,
         )
 
         mask = tailcuts_clean(

examples/camera_norms.py

@@ -22,16 +22,16 @@
 
     model = toymodel.generate_2d_shower_model(
         centroid=(0.2, 0.0),
-        width=0.01,
-        length=0.1,
+        width=0.05,
+        length=0.15,
         psi='35d',
     )
 
     image, sig, bg = toymodel.make_toymodel_shower_image(
         geom,
         model.pdf,
-        intensity=50,
-        nsb_level_pe=1000,
+        intensity=1500,
+        nsb_level_pe=5,
     )
 
     disps = []