From 3027265a1ba7e9df27859df222b82cb931f7d99a Mon Sep 17 00:00:00 2001 From: burmist-git Date: Tue, 24 Mar 2026 10:53:28 +0100 Subject: [PATCH 1/3] add polygon_chord --- src/ctapipe/image/muon/intensity_fitter.py | 36 ++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/src/ctapipe/image/muon/intensity_fitter.py b/src/ctapipe/image/muon/intensity_fitter.py index b3f3bb7e11f..ecadc3b387e 100644 --- a/src/ctapipe/image/muon/intensity_fitter.py +++ b/src/ctapipe/image/muon/intensity_fitter.py @@ -39,6 +39,42 @@ SQRT2 = np.sqrt(2) +def polygon_chord(mu_x, mu_y, phi, ri_x, ri_y, vi_x, vi_y): + vmu_x = np.cos(phi) + vmu_y = np.sin(phi) + + c1 = mu_x - ri_x + c2 = mu_y - ri_y + D = vi_x * vmu_y - vi_y * vmu_x + 1.0e-20 + + t = ( c1 * vmu_y - c2 * vmu_x ) / D + s = ( vi_y * c1 - vi_x * c2 ) / D + + status = np.column_stack((vi_x, ri_x, vi_y, ri_y, t, s)) + mask = (status[:,4] >= 0) & (status[:,4] < 1) & (status[:,5] >= 0) + + x_int = status[mask][:,0]*status[mask][:,4] + status[mask][:,1] + y_int = status[mask][:,2]*status[mask][:,4] + status[mask][:,3] + + if x_int.shape[0] == 0 : + return 0.0 + elif x_int.shape[0] == 1: + return np.squeeze(np.sqrt((x_int - mu_x) ** 2 + (y_int - mu_y) ** 2)) + elif x_int.shape[0] == 2: + return np.squeeze(np.sqrt((x_int[0] - x_int[1]) ** 2 + (y_int[0] - y_int[1]) ** 2)) + else: + dist = np.sort( + np.squeeze(np.sqrt((x_int - mu_x) ** 2 + (y_int - mu_y) ** 2)) + ) + sign_arr = np.ones(x_int.shape[0]) + if x_int.shape[0] % 2 == 0: + sign_arr[0::2] = -1 + return np.sum(dist * sign_arr) + + sign_arr[1::2] = -1 + return np.sum(dist * sign_arr) + + def chord_length(radius, rho, phi, phi0=0): """ Function for integrating the length of a chord across a circle (effective chord length). From 281155328fb41ad5bd1fb45cb0f4096297635700 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Mon, 30 Mar 2026 10:52:01 +0200 Subject: [PATCH 2/3] Add function documentation --- src/ctapipe/image/muon/intensity_fitter.py | 57 +++++++++++++++++++++- 1 file changed, 56 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/image/muon/intensity_fitter.py b/src/ctapipe/image/muon/intensity_fitter.py index ecadc3b387e..8c7a589f4aa 100644 --- a/src/ctapipe/image/muon/intensity_fitter.py +++ b/src/ctapipe/image/muon/intensity_fitter.py @@ -40,12 +40,67 @@ def polygon_chord(mu_x, mu_y, phi, ri_x, ri_y, vi_x, vi_y): + """ + Compute the chord length of a ray intersecting a polygon. + + This function calculates the length of the segment formed by the intersection + of a ray (defined by an origin and direction) with a polygon defined by its + edges. The polygon is represented parametrically using starting points and + direction vectors for each edge. + + Parameters + ---------- + mu_x : float + X-coordinate of the ray origin, which is the muon's impact point (x). + mu_y : float + Y-coordinate of the ray origin, which is the muon's impact point (y). + phi : float + Angle defining the direction of the ray. + ri_x : ndarray of shape (N,) + X-coordinates of the starting points of the polygon edges. + ri_y : ndarray of shape (N,) + Y-coordinates of the starting points of the polygon edges. + vi_x : ndarray of shape (N,) + X-components of the edge direction vectors. + vi_y : ndarray of shape (N,) + Y-components of the edge direction vectors. + + Returns + ------- + float + Length of the chord formed by the intersection of the ray with the polygon: + - 0.0 if there is no intersection, + - distance from the ray origin to the intersection point if only one intersection, + - distance between two intersection points if exactly two intersections, + - accumulated segment length for multiple intersections (handles complex polygons). + + Notes + ----- + - The ray is defined parametrically as: + (x, y) = (mu_x, mu_y) + s * (cos(phi), sin(phi)), with s >= 0 + - Each polygon edge is defined as: + (x, y) = (ri_x, ri_y) + t * (vi_x, vi_y), with 0 <= t < 1 + - The function computes intersections by solving a 2D linear system. + - A small epsilon (`1e-20`) is added to the denominator to avoid division by zero. + - For multiple intersections, distances are sorted and combined with alternating + signs to compute the total chord length (useful for non-convex polygons). + + Examples + -------- + >>> length = polygon_chord(mu_x=0, mu_y=0, phi=np.pi/4, + ... ri_x=..., ri_y=..., vi_x=..., vi_y=...) + >>> print(length) + """ + + # Effective speed of the ray, with unit norm. vmu_x = np.cos(phi) vmu_y = np.sin(phi) + epsilon_d = 1.0e-20 + c1 = mu_x - ri_x c2 = mu_y - ri_y - D = vi_x * vmu_y - vi_y * vmu_x + 1.0e-20 + D = vi_x * vmu_y - vi_y * vmu_x + epsilon_d t = ( c1 * vmu_y - c2 * vmu_x ) / D s = ( vi_y * c1 - vi_x * c2 ) / D From e3b0f46015ccead89d0fc70193826aa2acc834d5 Mon Sep 17 00:00:00 2001 From: burmist-git Date: Mon, 30 Mar 2026 11:54:24 +0200 Subject: [PATCH 3/3] rename D to determinant --- src/ctapipe/image/muon/intensity_fitter.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/image/muon/intensity_fitter.py b/src/ctapipe/image/muon/intensity_fitter.py index 8c7a589f4aa..f858b5bf27f 100644 --- a/src/ctapipe/image/muon/intensity_fitter.py +++ b/src/ctapipe/image/muon/intensity_fitter.py @@ -81,7 +81,7 @@ def polygon_chord(mu_x, mu_y, phi, ri_x, ri_y, vi_x, vi_y): - Each polygon edge is defined as: (x, y) = (ri_x, ri_y) + t * (vi_x, vi_y), with 0 <= t < 1 - The function computes intersections by solving a 2D linear system. - - A small epsilon (`1e-20`) is added to the denominator to avoid division by zero. + - A small epsilon_d (`1e-20`) is added to the denominator to avoid division by zero. - For multiple intersections, distances are sorted and combined with alternating signs to compute the total chord length (useful for non-convex polygons). @@ -90,6 +90,7 @@ def polygon_chord(mu_x, mu_y, phi, ri_x, ri_y, vi_x, vi_y): >>> length = polygon_chord(mu_x=0, mu_y=0, phi=np.pi/4, ... ri_x=..., ri_y=..., vi_x=..., vi_y=...) >>> print(length) + """ # Effective speed of the ray, with unit norm. @@ -100,10 +101,10 @@ def polygon_chord(mu_x, mu_y, phi, ri_x, ri_y, vi_x, vi_y): c1 = mu_x - ri_x c2 = mu_y - ri_y - D = vi_x * vmu_y - vi_y * vmu_x + epsilon_d + determinant = vi_x * vmu_y - vi_y * vmu_x + epsilon_d - t = ( c1 * vmu_y - c2 * vmu_x ) / D - s = ( vi_y * c1 - vi_x * c2 ) / D + t = ( c1 * vmu_y - c2 * vmu_x ) / determinant + s = ( vi_y * c1 - vi_x * c2 ) / determinant status = np.column_stack((vi_x, ri_x, vi_y, ri_y, t, s)) mask = (status[:,4] >= 0) & (status[:,4] < 1) & (status[:,5] >= 0)