Add get_pixel_id method #918

moritzhuetten · 2019-01-09T14:43:59Z

Added CameraGeometry.get_pixel_id() method to return the pixel number(s) corresponding to a (or an array of) x,y position(s), to be used like:

>>> geom = CameraGeometry.from_name("LSTCam")
>>> x, y = 0.1 * u.m, 0.1 * u.m
>>> geom.get_pixel_id(x, y)
138
>>> x, y = [0.1, 0.3] * u.m, [0.1, 0.4] * u.m
>>> geom.get_pixel_id(x, y)
array([138, 876])

maxnoe · 2019-01-09T14:51:32Z

Awesome.

How does this behave for points outside the camera?

maxnoe · 2019-01-09T14:52:04Z

ctapipe/instrument/camera.py

+        pixel_centers = np.stack([self.pix_x.to_value(u.m),self.pix_y.to_value(u.m)]).T
+        points_searched = np.dstack([x.to_value(u.m),y.to_value(u.m)])
+
+        kdtree = KDTree(pixel_centers)


I would suggest to build the kdtree in __init__, as it is used by more than one method.

agreed - right now we always make the assumption that you cannot modify the geometry (many functions cache a result). It could also go as a lazyproperty like the neighbor matrix, since most of the time you don't need it.

put it into lazyproperty in 7114978

maxnoe · 2019-01-09T14:52:58Z

ctapipe/instrument/camera.py

+        # pixel. If not, pos lies outside the camera. This approach does not need to know the  
+        # particular pixel shape.
+
+        borderpix_ids = np.where(self.get_border_pixel_mask())[0]


I think this can be vectorized by just using the mask

moritzhuetten · 2019-01-09T15:18:54Z

How does this behave for points outside the camera?

Returns -1 and prints a warning.

moritzhuetten · 2019-01-09T15:46:42Z

At the moment, the Outside-of-camera-check is exact (for arbitrary pixel shapes), but presumes that camera_pixels[0] lies in the center of the camera (which may not be always the case). A much simpler, but only approximate alternative would be to check if the position distance to the nearest pixel center is smaller than than the radius of the pixels' circumcircle.

moritzhuetten · 2019-01-09T17:01:01Z

Dropped assumption on central pixel in e95ff10, still being exact. However, the code might be much simpler and & readable if the inside-camera-check is only approximate. Or doing an exact check in a different, smarter way.

codecov · 2019-01-09T17:18:47Z

Codecov Report

Merging #918 into master will decrease coverage by 3.97%.
The diff coverage is 86.04%.

@@            Coverage Diff             @@
##           master     #918      +/-   ##
==========================================
- Coverage   82.31%   78.34%   -3.98%     
==========================================
  Files         186      191       +5     
  Lines       10716    10778      +62     
==========================================
- Hits         8821     8444     -377     
- Misses       1895     2334     +439

Impacted Files	Coverage Δ
ctapipe/instrument/tests/test_camera.py	`100% <100%> (ø)`	⬆️
ctapipe/instrument/camera.py	`93.17% <84.21%> (-2.64%)`	⬇️
ctapipe/io/tests/test_targetio_event_source.py	`6.5% <0%> (-93.5%)`	⬇️
ctapipe/io/tests/test_lsteventsource.py	`13.79% <0%> (-86.21%)`	⬇️
ctapipe/io/tests/test_sst1meventsource.py	`14.28% <0%> (-85.72%)`	⬇️
ctapipe/io/tests/test_nectarcameventsource.py	`16.66% <0%> (-83.34%)`	⬇️
ctapipe/io/targetioeventsource.py	`14.72% <0%> (-80.63%)`	⬇️
ctapipe/io/sources.py	`0% <0%> (-79.55%)`	⬇️
ctapipe/io/nectarcameventsource.py	`18.34% <0%> (-77.88%)`	⬇️
ctapipe/io/lsteventsource.py	`17.83% <0%> (-76.33%)`	⬇️
... and 46 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update c5f609b...de4ff00. Read the comment docs.

maxnoe · 2019-01-10T08:16:14Z

ctapipe/instrument/camera.py

+
+        """
+
+        pixel_centers = np.stack([self.pix_x.to_value(u.m), self.pix_y.to_value(u.m)]).T


np.column_stack gives directly what you want here without having to transpose.

done in e43e781

maxnoe · 2019-01-10T08:18:49Z

ctapipe/instrument/camera.py

+                   outside camera
+        '''
+
+        if not np.all(self.pix_area == self.pix_area[0], axis=0):


Use if np.any(self.pix_area != self.pix_area[0]). all always needs to test every value, any exits at the first True value it encounters.

Also this looks like it should be a @lazyproperty, like all_pixel_areas_equal or better name you can come up with.

also might want to use a tolerance, in case of any rounding error in case the area of the pixel is calculated somewhere. (or np.isclose() - is there an anyclose() function?),

np.any(~np.isclose(...

done in 704d7cc and 5c9c185. I marked the lazyproperty internal (leading underscore). Correct me if that's a wrong use-case

maxnoe · 2019-01-10T08:20:28Z

ctapipe/instrument/camera.py

+            for borderpix_id in borderpix_ids_in_list:
+                index = np.where(pixel_ids == borderpix_id)[0][0]
+                # compare with inside pixel:
+                xprime = points_searched[0][index, 0] \


for calculations, use implicit line continuations in parenthese instead of explicit line continuations, like this:

a = ( b + c )

Also, converting units is kind of expensive, so you should do self.pix_y[index].to_value(u.m) instead of converting all values and then only using one.

done in fd5d48f

maxnoe · 2019-01-10T08:25:31Z

ctapipe/instrument/camera.py

+        points_searched = np.dstack([x.to_value(u.m), y.to_value(u.m)])
+
+        kdtree = self.all_pixels_kdtree
+        dist, pixel_ids = kdtree.query(points_searched)


the kdtree..query has an option distance_upper_bound for where you could use the outer circle diameter of a pixel. this would immediately return n_pixel for points far away from any pixel.

At a certain size, it also might be faster to build a kdtree for the points that are looked for.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.cKDTree.query.html#scipy.spatial.cKDTree.query

distance_upper_bound implemented in a34e9bd. I left the precise check for now (only done for points between circumference and actual surface), but the corresponding code can simply be removed (is only called if points haven't been sorted out already by distance_upper_bound). For many queried points, distance_upper_bound seems a factor 3 faster than the current exact check.

At a certain size, it also might be faster to build a kdtree for the points that are looked for

Absolutely, but I fear this might turn into some over-engineering?

kosack · 2019-01-16T15:01:43Z

ctapipe/instrument/camera.py

@@ -518,6 +533,69 @@ def get_shower_coordinates(self, x, y, psi):

        return longi, trans

+    def get_pixel_id(self, x, y):


I'd suggest renaming this - it's not clear what the input is ("get" could mean just lookup the pixel_id from the pixel index). Maybe position_to_pixel_id(x,y)?

renamed method to position_to_pix_index() in 49e3572

kosack · 2019-01-16T15:03:21Z

ctapipe/instrument/camera.py

+
+        Returns
+        -------
+        pixel_ids: Pixel number or array of pixel numbers. Returns -1 if position falls 


Does this give the pixel ID (the name of the pixel), or the pixel index the index into the arrays in the CameraGeometry (these are not the same thing)

given a pixel index i, the x-position of that pixel isgeom.pix_x[i], and the pixel_id of that pixel is geom.pix_id[i].

You should avoid the word "pixel number" and don't confuse the concepts.

Looking at the code, it seem to compute the pixel_index, not the pixel_id

Indeed, it's the pixel index. Wasn't aware of the difference in ctapipe, but of course makes sense to differentiate, thanks for pointing this out. Renamed the method to position_to_pix_index(), because thought it's useful to let the user decide if s/he wants to obtain the id in a separate step. Also clarified the variable naming to avoid confusing expression 'number'.

kosack · 2019-01-16T15:05:37Z

ctapipe/instrument/camera.py

+                   outside camera
+        '''
+
+        if not np.all(self.pix_area == self.pix_area[0], axis=0):


also might want to use a tolerance, in case of any rounding error in case the area of the pixel is calculated somewhere. (or np.isclose() - is there an anyclose() function?),

kosack · 2019-01-16T15:06:29Z

ctapipe/instrument/camera.py

@@ -203,6 +203,21 @@ def _calc_pixel_area(pix_x, pix_y, pix_type):

        return np.ones(pix_x.shape) * area

+    @lazyproperty
+    def all_pixels_kdtree(self):


is there a case where this is not for all pixels? Otherwise, perhaps just name this property kdtree

Renamed it to simply kdtree in f1a7b92. If finally satisfied with this PR, maybe have a look at already existing _find_neighbor_pixels method() to work with the same kdtree (seems to also use all pixels, but slightly differently implemented)

kosack · 2019-01-16T15:56:48Z

Dropped assumption on central pixel in e95ff10, still being exact. However, the code might be much simpler and & readable if the inside-camera-check is only approximate. Or doing an exact check in a different, smarter way.

Couldn't that just be a bounding box? E.g. (pix_x.min()-pix_rad, pix_x.max()+pix_rad), that's an easy and rough way to avoid computing if it's in the camera. You could compute a bounding circle as well, which might be better since most cameras are pretty circular, so r < hypot(pix_x, pix_y).max() or similar.

moritzhuetten · 2019-01-25T00:17:53Z

Dropped assumption on central pixel in e95ff10, still being exact. However, the code might be much simpler and & readable if the inside-camera-check is only approximate. Or doing an exact check in a different, smarter way.

Couldn't that just be a bounding box? E.g. (pix_x.min()-pix_rad, pix_x.max()+pix_rad), that's an easy and rough way to avoid computing if it's in the camera. You could compute a bounding circle as well, which might be better since most cameras are pretty circular, so r < hypot(pix_x, pix_y).max() or similar.

I implemented @maxnoe's suggestion to do it approximately with the circumference on the individual pixel basis. As commented above, can easily simply remove the complicated chunk of code now.

moritzhuetten added 2 commits January 9, 2019 13:47

camera method to return pixel of (x,y) position

1725800

cleaned up new get_pixel_id method + test,comments added

f7baf78

maxnoe reviewed Jan 9, 2019

View reviewed changes

moritzhuetten added 2 commits January 9, 2019 17:51

put all camera pixels' kdtree into lazyproperty function

7114978

drop assumption camera_pixels[0] not at the border

e95ff10

Comparison pixel selection inside condition

70419c4

moritzhuetten added 2 commits January 10, 2019 08:15

warning added that kdtree approach only works for equal pixel sizes

ff6519a

PEP8 compliant code

c33465a

maxnoe reviewed Jan 10, 2019

View reviewed changes

kosack requested changes Jan 16, 2019

View reviewed changes

moritzhuetten added 9 commits January 24, 2019 20:03

use numpy.column_stack() instead of transpose

e43e781

use numpy.any() and np.isclose() to check equal pixel sizes

704d7cc

renamed method to position_to_pix_index()

49e3572

renamed all_pixels_kdtree to kdtree

f1a7b92

created lazyproperty _all_pixel_areas_equal()

5c9c185

updated test_camera.py file

36f06b5

implemented circumference-limited kdtree + code style improvement

a34e9bd

only convert units for single array entry

fd5d48f

removed trailing whitespace

5215c72

tests should pass now

de4ff00

kosack approved these changes Jan 28, 2019

View reviewed changes

dneise approved these changes Jan 28, 2019

View reviewed changes

dneise merged commit 6bc409b into cta-observatory:master Jan 28, 2019

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add get_pixel_id method #918

Add get_pixel_id method #918

moritzhuetten commented Jan 9, 2019

maxnoe commented Jan 9, 2019

maxnoe Jan 9, 2019

kosack Jan 9, 2019

moritzhuetten Jan 9, 2019

maxnoe Jan 9, 2019

moritzhuetten commented Jan 9, 2019

moritzhuetten commented Jan 9, 2019

moritzhuetten commented Jan 9, 2019

codecov bot commented Jan 9, 2019 •

edited

Loading

maxnoe Jan 10, 2019

moritzhuetten Jan 24, 2019 •

edited

Loading

maxnoe Jan 10, 2019

kosack Jan 16, 2019

maxnoe Jan 16, 2019

moritzhuetten Jan 24, 2019 •

edited

Loading

maxnoe Jan 10, 2019

maxnoe Jan 10, 2019

moritzhuetten Jan 25, 2019 •

edited

Loading

maxnoe Jan 10, 2019 •

edited

Loading

maxnoe Jan 10, 2019

moritzhuetten Jan 24, 2019 •

edited

Loading

kosack Jan 16, 2019

moritzhuetten Jan 24, 2019 •

edited

Loading

kosack Jan 16, 2019

kosack Jan 16, 2019

moritzhuetten Jan 24, 2019

kosack Jan 16, 2019

kosack Jan 16, 2019

moritzhuetten Jan 24, 2019 •

edited

Loading

kosack commented Jan 16, 2019

moritzhuetten commented Jan 25, 2019


		"""

		pixel_centers = np.stack([self.pix_x.to_value(u.m), self.pix_y.to_value(u.m)]).T

		@@ -518,6 +533,69 @@ def get_shower_coordinates(self, x, y, psi):

		return longi, trans

		def get_pixel_id(self, x, y):

Add get_pixel_id method #918

Add get_pixel_id method #918

Conversation

moritzhuetten commented Jan 9, 2019

maxnoe commented Jan 9, 2019

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

moritzhuetten commented Jan 9, 2019

moritzhuetten commented Jan 9, 2019

moritzhuetten commented Jan 9, 2019

codecov bot commented Jan 9, 2019 • edited Loading

Codecov Report

Choose a reason for hiding this comment

moritzhuetten Jan 24, 2019 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

moritzhuetten Jan 24, 2019 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

moritzhuetten Jan 25, 2019 • edited Loading

Choose a reason for hiding this comment

maxnoe Jan 10, 2019 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

moritzhuetten Jan 24, 2019 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

moritzhuetten Jan 24, 2019 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

moritzhuetten Jan 24, 2019 • edited Loading

Choose a reason for hiding this comment

kosack commented Jan 16, 2019

moritzhuetten commented Jan 25, 2019

codecov bot commented Jan 9, 2019 •

edited

Loading

moritzhuetten Jan 24, 2019 •

edited

Loading

moritzhuetten Jan 24, 2019 •

edited

Loading

moritzhuetten Jan 25, 2019 •

edited

Loading

maxnoe Jan 10, 2019 •

edited

Loading

moritzhuetten Jan 24, 2019 •

edited

Loading

moritzhuetten Jan 24, 2019 •

edited

Loading

moritzhuetten Jan 24, 2019 •

edited

Loading