add polygon_chord by burmist-git · Pull Request #2969 · cta-observatory/ctapipe

burmist-git · 2026-03-24T09:55:35Z

Polygonal description of the chord.

096_arbitrary_shaped_mirror_and_polygon_chord.pdf

Analytical solution :

Input:
Ray origin, ray direction, polygon edges

Compute the ray direction vector.
Compute the intersection parameters (t, s) for all polygon edges.
Keep only intersections satisfying:
0 ≤ t < 1 and s ≥ 0.
Compute the intersection coordinates.
If no intersections:
return 0
If one intersection:
return distance from the ray origin to the intersection
If two intersections:
return distance between the two intersections
Otherwise:
Sort the distances from the ray origin to all intersections.
Return the alternating sum of the sorted distances.

Approximation with Fourier transformation:

kosack · 2026-03-24T17:57:13Z

Maybe we should consider adding shapely as a dependency? It might simplify some of these computations without having to re-implement all of computational geometry.

maxnoe · 2026-03-26T15:09:21Z

Maybe we should consider adding shapely as a dependency?

I already looked into using shapely for a way to support trigger / voltage patches in the past here:
#855 (comment)

burmist-git · 2026-03-26T18:53:33Z

This is possible, but the function itself is very small, so using an external library would be unnecessary. I wouldn’t make it more complex than it is.

mexanick · 2026-03-31T07:11:40Z

+        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,)


I think we should deduce these (if needed) from the vertices. You can write a helper for this.

Yes, I'd also prefer the vertices as input. But maybe the function as it is here is fine and that transformation is made once outside of this function.

Yes agree with @maxnoe . This is a low-level function. Let's keep it as it is, since all the necessary computations should be handled by the higher-level function that calls it.

mexanick · 2026-03-31T07:14:32Z

+
+    Parameters
+    ----------
+    mu_x : float


These should be quantities, as we discussed offline, the best would be to use telescope coordinate system lat/lon.

Also, perhaps tuples of (x,y) or arrays of them would be more straightforward to define points (arrays of points)

mexanick · 2026-03-31T07:15:20Z

+
+    Returns
+    -------
+    float


Should be also a unit of length (or angular, equivalent in the telescope frame)

mexanick · 2026-03-31T07:19:42Z

+    - 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_d (`1e-20`) is added to the denominator to avoid division by zero.


why not using

try: ... except ZeroDivisionError: # determinant is zero, the cord goes on the edge, handle this

Just remember that catching an exception is quite a bit slower than adding an epsilon, so if you need this to be fast, the epsilon method might be better. Also if you want to use numba to speed things up

And doesn't vectorize

maxnoe · 2026-03-31T07:50:57Z

 SQRT2 = np.sqrt(2)


+def polygon_chord(mu_x, mu_y, phi, ri_x, ri_y, vi_x, vi_y):


This function should be njit, as it will be called. If possible, it would also be good if it could be vectorized over multiple tiles.

Yes, I see. Let me finish implementing it with multiple tiles.

Voutsi · 2026-08-03T06:13:00Z

+
+    Returns
+    -------
+    float


Every branch returns (length, x_int, y_int) and not only a float. This should be reflected in the docstring.

It is indeed only a float. The function is not vectorized.

Voutsi · 2026-08-03T06:18:35Z

+    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 :


I suggest to return the same type (an array) for every if loop branch. I don't say that the code is wrong, but returning the same type will make it cleaner.

It is indeed only a float. The function is not vectorized. All of them return only one number.

… set of projection angles.

ctao-sonarqube · 2026-08-07T09:33:15Z

Quality Gate failed

Failed conditions
1 New issue
74.1% Coverage on New Code (required ≥ 80%)

See analysis details on SonarQube

Catch issues before they fail your Quality Gate with our IDE extension SonarQube for IDE

burmist-git · 2026-08-10T16:15:59Z

@Voutsi @kosack @maxnoe The code, tests, and documentation are ready. All tests are green and have passed. There is currently an issue with the SonarQube check related to code complexity (there are several nested if conditions). I know how to fix it, but before making the change, I want to make sure that the current architecture is appropriate and will be accepted for merging. I will fix it as soon as I receive feedback from the reviewers. For the tests with the full simulation chain (not the unit tests), I have requested the simulation from SimPipe.https://gitlab.cta-observatory.org/groups/cta-computing/dpps/-/work_items?sort=created_date&state=opened&search=muons&first_page_size=20&show=eyJpaWQiOiIyOTUiLCJmdWxsX3BhdGgiOiJjdGEtY29tcHV0aW5nL2RwcHMvY2FsaWJyYXRpb25waXBlbGluZS9jYWxpYnBpcGUiLCJpZCI6MTg4OTJ9

maxnoe · 2026-08-10T16:35:41Z

The code, tests, and documentation are ready.

The polygon_chord function is a low-level function. I'd expect it to be (maybe optionally) used in the muon fitter.

At the moment, this function is not used anywhere. I don't really see a benefit in adding this function without it being possible to actually use it in the muon fit.

maxnoe · 2026-08-10T16:36:43Z

I also thought that we had concluded at the meeting in Geneva that we would need to handle individual mirror tiles, so what is the motivation for still adding this polygon approximation?

maxnoe · 2026-08-10T16:38:32Z

+    vi_y = []
+
+    for ver_i in vertices_list:
+        ver_f = np.roll(ver_i, 1, axis=0)


In my experience, np.roll is extremely expensive, especially when called in a loop here, as it makes full copies of its input data.

In general I am much worried about the (runtime) performance of the code here, as it contains a lot of raw python loops without numba compilation for functions that are meant to be called in a likelihood optimization.

For the optimization, let’s leave it for the next iteration. The Numba optimization can be applied to cases with zero, one, or two intersection points. However, for more complex shapes, such as the one we want to have for a the LST mirror (not the single facet), this requires a sorting function that cannot be used with Numba.

maxnoe · 2026-08-10T16:40:10Z

+                np.squeeze(x_int),
+                np.squeeze(y_int),
+            )
+        return np.squeeze(np.sqrt((x_int - mu_x) ** 2 + (y_int - mu_y) ** 2))


please do not duplicate computation of the same things in different branches.

compute first, then return, i.e.

chord_length = ... if return intersections: return chord_length, ... return chord_length

burmist-git · 2026-08-10T16:40:51Z

I also thought that we had concluded at the meeting in Geneva that we would need to handle individual mirror tiles, so what is the motivation for still adding this polygon approximation?

The hexagonal shape of the individual tiles can also be used or any other shape. However, that was not really the conclusion we reached regarding the use of the tiles individually. We considered using the individual tiles as an option, and this is now possible.

burmist-git · 2026-08-10T16:43:04Z

The code, tests, and documentation are ready.

The polygon_chord function is a low-level function. I'd expect it to be (maybe optionally) used in the muon fitter.

At the moment, this function is not used anywhere. I don't really see a benefit in adding this function without it being possible to actually use it in the muon fit.

Sure, this has to be done. I would consider adding it in the next PR, but if needed, I can add it now.

maxnoe · 2026-08-10T16:44:41Z

As I am very worried about the performance of the code being added here (I think it will be prohibitively slow as is), I'd like to see the performance of it when used as part of the fitter before merging.

burmist-git · 2026-08-10T16:51:35Z

As I am very worried about the performance of the code being added here (I think it will be prohibitively slow as is), I'd like to see the performance of it when used as part of the fitter before merging.

These are the requirements - we need to fulfill them. As we agreed during the meeting, performance comes later as a requirement. So, this option should be available in ctapipe.

maxnoe · 2026-08-10T16:54:34Z

I disagree, if we see here that the performance is as bad as I expect it to be from looking at the code here, this will not be a question of CalibPipe performance requirements, but just that I won't merge this in ctapipe.

burmist-git · 2026-08-10T16:57:52Z

I disagree, if we see here that the performance is as bad as I expect it to be from looking at the code here, this will not be a question of CalibPipe performance requirements, but just that I won't merge this in ctapipe.

Please add this comment here, under the Level C requirements for muons : https://gitlab.cta-observatory.org/cta-computing/dpps/calibrationpipeline/documentation/calibpipe-level-c-requirements/-/merge_requests/10

| C-CALP-0308 | Arbitrary description of the primary mirror shape. | CalibPipe shall account for the non-circularity of the reflectors when reconstructing the muon impact point. | T: Dedicated simulations of muon events with a fixed impact point shall be processed using CalibPipe. Qualification shall be conducted by comparing the reconstructed impact point against the ground truth, taking into account the corresponding uncertainties. | CTAO-B_DPPS-494, CTAO-B_DPPS-439, CTAO-B_DPPS-408 | UC-120.2.2, UC-120.2.2.1, UC-120.2.2.2, UC-120.2.2.3, UC-120.2.10, UC-120.2.10.1, UC-120.2.10.2, UC-120.2.10.3, UC-120.2.13, UC-120.2.14|

Now we are in a position to keep track of whether the requirements have been fulfilled or not, and why.

burmist-git · 2026-08-10T19:36:47Z

@maxnoe please look at the functional requirements regarding C-CALP-0308 : https://gitlab.cta-observatory.org/cta-computing/dpps/calibrationpipeline/documentation/calibpipe-level-c-requirements/-/merge_requests/10#note_422640

burmist-git · 2026-08-11T08:45:04Z

I understand the overall feedback and will work on optimizing the code and incorporating it into the intensity fitter.

Voutsi · 2026-08-13T08:14:40Z

I understand the overall feedback and will work on optimizing the code and incorporating it into the intensity fitter.

Thanks @burmist-git

I disagree, if we see here that the performance is as bad as I expect it to be from looking at the code here

@maxnoe , will be the optimisation discussed above enough in order for this PR to be acceptable? I am speaking only performance-wise, of course all other threads should be as well resolved.

burmist-git · 2026-08-13T08:36:52Z

Just in case, @maxnoe, this is the optimization I’m proposing :

#2969 (comment)

maxnoe · 2026-08-13T08:52:28Z

+    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))


Why do you column_stack this (making a copy) just to then index into the individual comments making it basically unreadable what you are accessing?

Yes, thank you. It can be removed. Initially, I did it to have everything in the same 2D array, but you’re right, it can be removed.

add polygon_chord

3027265

burmist-git requested a review from maxnoe March 24, 2026 09:56

burmist-git added 2 commits March 30, 2026 10:52

Add function documentation

2811553

rename D to determinant

e3b0f46

mexanick reviewed Mar 31, 2026

View reviewed changes

maxnoe reviewed Mar 31, 2026

View reviewed changes

burmist-git mentioned this pull request Apr 13, 2026

Non-spherical reflectors for muon analysis. #2988

Open

burmist-git requested a review from Voutsi July 13, 2026 14:02

burmist-git self-assigned this Jul 13, 2026

add the intersection points

a64959b

Voutsi reviewed Aug 3, 2026

View reviewed changes

burmist-git added 15 commits August 4, 2026 11:30

Add debug flag return_intersections=False

43efa97

Add debug flag return_intersections=False

9987b67

add Compute the total chord length through one or more polygons for a…

5d099f0

… set of projection angles.

Add test_polygon_chord

31ba178

Add change log

8b8f76a

add polygon_chord_base

f067c07

add polygon_chord_base

1303aae

Add test_polygon_chord_base and test_polygon_chord

2805134

Update

788c296

fix style

3919025

fix style

5995c7e

fix style

b287196

fix style

b1f2387

fix style

ffeb580

fix style

b7c28ae

burmist-git added 4 commits August 5, 2026 17:33

fix style

36ad270

ruff format intensity_fitter.py

ca8802d

Update

3c48d60

Update doc string

d46c1a4

burmist-git requested review from Voutsi and kosack August 5, 2026 16:23

fix doc string

0b92660

burmist-git requested a review from maxnoe August 10, 2026 08:19

maxnoe reviewed Aug 10, 2026

View reviewed changes

maxnoe reviewed Aug 13, 2026

View reviewed changes

		SQRT2 = np.sqrt(2)


		def polygon_chord(mu_x, mu_y, phi, ri_x, ri_y, vi_x, vi_y):

Uh oh!

Conversation

burmist-git commented Mar 24, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

kosack commented Mar 24, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

maxnoe commented Mar 26, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

burmist-git commented Mar 26, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

kosack Apr 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ctao-sonarqube Bot commented Aug 7, 2026

Quality Gate failed

Uh oh!

burmist-git commented Aug 10, 2026

Uh oh!

maxnoe commented Aug 10, 2026

Uh oh!

maxnoe commented Aug 10, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

burmist-git commented Aug 10, 2026

Uh oh!

burmist-git commented Aug 10, 2026

Uh oh!

maxnoe commented Aug 10, 2026

Uh oh!

burmist-git commented Aug 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

maxnoe commented Aug 10, 2026

Uh oh!

burmist-git commented Aug 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

burmist-git commented Aug 10, 2026

Uh oh!

burmist-git commented Aug 11, 2026

Uh oh!

Voutsi commented Aug 13, 2026

Uh oh!

burmist-git commented Mar 24, 2026 •

edited

Loading

kosack commented Mar 24, 2026 •

edited

Loading

maxnoe commented Mar 26, 2026 •

edited

Loading

kosack Apr 13, 2026 •

edited

Loading

burmist-git commented Aug 10, 2026 •

edited

Loading

burmist-git commented Aug 10, 2026 •

edited

Loading