Fix hillas width 0, fixes #772

ctapipe/image/hillas.py

@@ -11,6 +11,9 @@
 from ..io.containers import HillasParametersContainer
 
 
+HILLAS_ATOL = np.finfo(np.float64).eps
+
+
 __all__ = [
     'hillas_parameters',
     'HillasParameterizationError',
@@ -137,6 +140,10 @@ def hillas_parameters(geom, image):
     cov = np.cov(delta_x, delta_y, aweights=image, ddof=0)
     eig_vals, eig_vecs = np.linalg.eigh(cov)
 
+    # round eig_vals to get rid of nans when eig val is something like -8.47032947e-22
+    near_zero = np.isclose(eig_vals, 0, atol=HILLAS_ATOL)
+    eig_vals[near_zero] = 0
+
     # width and length are eigen values of the PCA
     width, length = np.sqrt(eig_vals)
 

ctapipe/image/tests/test_hillas.py

@@ -202,3 +202,35 @@ def test_skewness():
             assert result.skewness == approx(-skew, abs=0.3)
 
         assert signal.sum() == result.intensity
+
+
+def test_straight_line_width_0():
+    ''' Test that hillas_parameters.width is 0 for a straight line of pixels '''
+    # three pixels in a straight line
+    long = np.array([0, 1, 2]) * 0.01
+    trans = np.zeros(len(long))
+    pix_id = np.arange(len(long))
+
+    np.random.seed(0)
+
+    for dx in (-1, 0, 1):
+        for dy in (-1, 0, 1):
+            for psi in np.linspace(0, np.pi, 20):
+                x = dx + np.cos(psi) * long + np.sin(psi) * trans
+                y = dy - np.sin(psi) * long + np.cos(psi) * trans
+
+                geom = CameraGeometry(
+                    cam_id='testcam',
+                    pix_id=pix_id,
+                    pix_x=x * u.m,
+                    pix_y=y * u.m,
+                    pix_type='hexagonal',
+                    pix_area=1 * u.m**2,
+                    sampling_rate=u.Quantity(1, u.GHz),
+                    reference_pulse_shape=np.ones(1),
+                    reference_pulse_step=u.Quantity(1, u.ns),
+                )
+
+                img = np.random.poisson(5, size=len(long))
+                result = hillas_parameters(geom, img)
+                assert result.width.value == 0

docs/tutorials/theta_square.ipynb

@@ -132,9 +132,13 @@
     "        # Calulate hillas parameters\n",
     "        # It fails for empty pixels\n",
     "        try:\n",
-    "            hillas_params[tel_id] = hillas_parameters(camgeom, cleaned_image)\n",
-    "        except:\n",
+    "            hillas = hillas_parameters(camgeom, cleaned_image)\n",
+    "        except Exception as e:\n",
+    "            print(e)\n",
     "            pass\n",
+    "        \n",
+    "        if hillas.width.value > 0:\n",
+    "            hillas_params[tel_id] = hillas\n",
     "\n",
     "    if len(hillas_params) < 2:\n",
     "        continue\n",
@@ -217,13 +221,6 @@
    "source": [
     "again, this plot is not beautiful since we have such low stats"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -242,7 +239,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.2"
+   "version": "3.8.1"
   },
   "toc": {
    "nav_menu": {