Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Drop support for CameraFrame in ImageProcessor #2141

Closed
wants to merge 23 commits into from
Closed
Show file tree
Hide file tree
Changes from 20 commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
d0b8083
remove CameraFrame from hillas
StFroese Nov 28, 2022
8034cee
remove CameraHillasParametersContainer option
StFroese Nov 28, 2022
a004399
extraction in TelescopeFrame
StFroese Nov 28, 2022
430fed2
remove CameraHillasParametersContainer option
StFroese Nov 28, 2022
26068fd
toymodels in TelescopeFrame
StFroese Nov 28, 2022
ed0f266
tests in TelescopeFrame
StFroese Nov 28, 2022
a0a9ee3
image processing in TelescopeFrame only
StFroese Nov 28, 2022
61005fd
remove CameraFrame option
StFroese Nov 28, 2022
67a9143
test in TelescopeFrame
StFroese Nov 28, 2022
765e50b
remove test for CameraFrame
StFroese Nov 28, 2022
a6eb055
create ToyEventSource in TelescopeFrame
StFroese Nov 28, 2022
3440e85
remove comparison between CameraFrame and TelescopeFrame
StFroese Nov 28, 2022
9b93b91
remove CameraFrame
StFroese Nov 28, 2022
37fe0e7
remove CameraHillasParametersContainer
StFroese Dec 7, 2022
2bc2b3e
speedup using frame instead of SkyCoord
StFroese Dec 8, 2022
9af1cf7
Merge branch 'master' into remove_hillas_cameraframe
StFroese Nov 23, 2022
2f3c6de
pixel_positions for TelescopeFrame
StFroese Jan 9, 2023
712e436
use TelescopeFrame in doc notebooks
StFroese Jan 9, 2023
84ceeb7
Merge branch 'master' into remove_hillas_cameraframe
StFroese Jan 9, 2023
de32018
fix CameraDisplay test
StFroese Jan 11, 2023
72c33b2
move geometry transformation to setup
StFroese Jan 30, 2023
3128997
subarray transform to TelescopeFrame
StFroese Jan 30, 2023
ce99544
Merge branch 'master' into remove_hillas_cameraframe
StFroese Jan 30, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 5 additions & 35 deletions ctapipe/containers.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,8 +32,6 @@
"MonitoringCameraContainer",
"MonitoringContainer",
"MorphologyContainer",
"BaseHillasParametersContainer",
"CameraHillasParametersContainer",
"CameraTimingParametersContainer",
"ParticleClassificationContainer",
"PedestalContainer",
Expand Down Expand Up @@ -184,45 +182,17 @@ class TelEventIndexContainer(Container):
tel_id = tel_id_field()


class BaseHillasParametersContainer(Container):
"""
Base container for hillas parameters to
allow the CameraHillasParametersContainer to
be assigned to an ImageParametersContainer as well.
"""

intensity = Field(nan, "total intensity (size)")
skewness = Field(nan, "measure of the asymmetry")
kurtosis = Field(nan, "measure of the tailedness")


class CameraHillasParametersContainer(BaseHillasParametersContainer):
"""
Hillas Parameters in the camera frame. The cog position
is given in meter from the camera center.
"""

default_prefix = "camera_frame_hillas"
x = Field(nan * u.m, "centroid x coordinate", unit=u.m)
y = Field(nan * u.m, "centroid x coordinate", unit=u.m)
r = Field(nan * u.m, "radial coordinate of centroid", unit=u.m)
phi = Field(nan * u.deg, "polar coordinate of centroid", unit=u.deg)

length = Field(nan * u.m, "standard deviation along the major-axis", unit=u.m)
length_uncertainty = Field(nan * u.m, "uncertainty of length", unit=u.m)
width = Field(nan * u.m, "standard spread along the minor-axis", unit=u.m)
width_uncertainty = Field(nan * u.m, "uncertainty of width", unit=u.m)
psi = Field(nan * u.deg, "rotation angle of ellipse", unit=u.deg)


class HillasParametersContainer(BaseHillasParametersContainer):
class HillasParametersContainer(Container):
"""
Hillas Parameters in a spherical system centered on the pointing position
(TelescopeFrame). The cog position is given as offset in
longitude and latitude in degree.
"""

default_prefix = "hillas"
intensity = Field(nan, "total intensity (size)")
skewness = Field(nan, "measure of the asymmetry")
kurtosis = Field(nan, "measure of the tailedness")
fov_lon = Field(
nan * u.deg,
"longitude angle in a spherical system centered on the pointing position",
Expand Down Expand Up @@ -366,7 +336,7 @@ class ImageParametersContainer(Container):
hillas = Field(
default_factory=HillasParametersContainer,
description="Hillas Parameters",
type=BaseHillasParametersContainer,
type=HillasParametersContainer,
)
timing = Field(
default_factory=TimingParametersContainer,
Expand Down
50 changes: 17 additions & 33 deletions ctapipe/image/concentration.py
Original file line number Diff line number Diff line change
@@ -1,14 +1,10 @@
import numpy as np
import astropy.units as u
import numpy as np

from ..containers import (
ConcentrationContainer,
HillasParametersContainer,
CameraHillasParametersContainer,
)
from .hillas import camera_to_shower_coordinates
from ..containers import ConcentrationContainer
from ..instrument import CameraGeometry
from ..utils.quantities import all_to_value
from .hillas import camera_to_shower_coordinates

__all__ = ["concentration_parameters"]

Expand All @@ -24,44 +20,32 @@ def concentration_parameters(geom: CameraGeometry, image, hillas_parameters):
"""

h = hillas_parameters
if isinstance(h, CameraHillasParametersContainer):
unit = h.x.unit
pix_x, pix_y, x, y, length, width, pixel_width = all_to_value(
geom.pix_x,
geom.pix_y,
h.x,
h.y,
h.length,
h.width,
geom.pixel_width,
unit=unit,
)
elif isinstance(h, HillasParametersContainer):
unit = h.fov_lon.unit
pix_x, pix_y, x, y, length, width, pixel_width = all_to_value(
geom.pix_x,
geom.pix_y,
h.fov_lon,
h.fov_lat,
h.length,
h.width,
geom.pixel_width,
unit=unit,
)

unit = h.fov_lon.unit
pix_x, pix_y, x, y, length, width, pixel_width = all_to_value(
geom.pix_x,
geom.pix_y,
h.fov_lon,
h.fov_lat,
h.length,
h.width,
geom.pixel_width,
unit=unit,
)

delta_x = pix_x - x
delta_y = pix_y - y

# take pixels within one pixel diameter from the cog
mask_cog = (delta_x ** 2 + delta_y ** 2) < pixel_width ** 2
mask_cog = (delta_x**2 + delta_y**2) < pixel_width**2
conc_cog = np.sum(image[mask_cog]) / h.intensity

if hillas_parameters.width.value != 0:
# get all pixels inside the hillas ellipse
longi, trans = camera_to_shower_coordinates(
pix_x, pix_y, x, y, h.psi.to_value(u.rad)
)
mask_core = (longi ** 2 / length ** 2) + (trans ** 2 / width ** 2) <= 1.0
mask_core = (longi**2 / length**2) + (trans**2 / width**2) <= 1.0
conc_core = image[mask_core].sum() / h.intensity
else:
conc_core = 0.0
Expand Down
11 changes: 9 additions & 2 deletions ctapipe/image/extractor.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
from traitlets import Bool, Int

from ctapipe.containers import DL1CameraContainer
from ctapipe.coordinates import TelescopeFrame
from ctapipe.core import TelescopeComponent
from ctapipe.core.traits import (
BoolTelescopeParameter,
Expand Down Expand Up @@ -1051,7 +1052,9 @@ def _apply_second_pass(
# Apply correction to 1st pass charges
charge_1stpass = charge_1stpass_uncorrected * correction[selected_gain_channel]

camera_geometry = self.subarray.tel[tel_id].camera.geometry
camera_geometry = self.subarray.tel[tel_id].camera.geometry.transform_to(
TelescopeFrame()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this might need some thought - we shouldn't be transforming the geometry for each event inside a calibration/image-extraction algorithm like this as that adds quite a bit of overhead. This algorithm does not care about the frame of the camera geometry or pointing corrections, so it should not require a transform.

Since we do eventually need to transform for the ImageProcessor, it may be better to do the transform outside of these algorithms, perhaps in the CameraCalibrator (which eventually should also apply pointing correcions), so that a single TelescopeFrame is available to all algorithms that need it. That would avoid re-doing the transform many times.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this should be done inside the ProcessorTool in setup after writing out the subarray description but before setting up telescope components. The telescope components should just get a subarray with cameras in CameraFrame.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe add transform_camera_geometries to SubarrayDescription returning a new Subarray with transformed geometries

Copy link
Contributor

@kosack kosack Jan 16, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is creating a whole new SubarrayDescription each time efficient? This was an issue I think we discussed a while back (need to find where): during analysis, the TelescopeFrame needs to be updated once per event for observations, and could be simplified to only once per run for simulations in fixed Alt/Az mode.

The TelescopeFrame requires:

  • the latest pointing position, which changes in time in ICRS tracking mode (the default for most observations)
  • pre-computed pointing corrections which could include not just a shift in RA/Dec, but also a rotation (this is not currently implemented, but will need to be done eventually)

Therefore where to store these frames is a question. Normally the SubarrayDescription is a static thing, but it could go there, but re-generating a whole SubarrayDescription might be wasteful if we could just update the info. Note that we don't want to re-generate neighbor pixel info in CameraGeometry for each event, as that does not change.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is creating a whole new SubarrayDescription each time efficient?

Not each time. Once in the tool in setup.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, sorry, forgot about pointing corrections.

But it stands: we could also do the pointing corrections in the transformation from / to TelescopeFrame to higher frames. That would actually make the most sense to me, since pointing direction is not important for the CameraFrame <-> TelescopeFrame transformation.

To me that means we should transform the whole subarray once to TelescopeFrame and apply corrections in the transformation AltAz <-> TelescopeFrame, we don't need to change telescope frame for this, as the attributes needed for the transformation can be specified on the SkyCoord as seen here:

coord : `~astropy.coordinates.SkyCoord`
The coordinate to plot. Must be able to be transformed into
the frame of the camera geometry of this display.
Most of the time, this means you need to add the telescope
pointing as a coordinate attribute like this:
``SkyCoord(..., telescope_pointing=pointing)``

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remember there are in principle rotation corrections and maybe others, so you have to apply pointing corrections before Hillas parameterization (doing it while going to Alt/Az wouldn't work), as you need to have correctly shifted/rotated parameters before you do shower geometry reconstruction. But that can be solved when the ShowerProcessor is moved to TelescopeFrame later (right now e.g. HillasReconstructor supports both frames, but that's ugly).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@maxnoe @kosack please have a look now if this solves the issue

)
if self.invalid_pixel_handler is not None:
charge_1stpass, pulse_time_1stpass = self.invalid_pixel_handler(
tel_id,
Expand Down Expand Up @@ -1122,7 +1125,11 @@ def _apply_second_pass(

# get projected distances along main image axis
longitude, _ = camera_to_shower_coordinates(
camera_geometry.pix_x, camera_geometry.pix_y, hillas.x, hillas.y, hillas.psi
camera_geometry.pix_x,
camera_geometry.pix_y,
hillas.fov_lon,
hillas.fov_lat,
hillas.psi,
)

# get the predicted times as a linear relation
Expand Down
19 changes: 2 additions & 17 deletions ctapipe/image/hillas.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
import numpy as np
from astropy.coordinates import Angle

from ..containers import CameraHillasParametersContainer, HillasParametersContainer
from ..containers import HillasParametersContainer

HILLAS_ATOL = np.finfo(np.float64).eps

Expand Down Expand Up @@ -73,7 +73,7 @@ def hillas_parameters(geom, image):

Parameters
----------
geom: ctapipe.instrument.CameraGeometry
geom: ctapipe.instrument.CameraGeometry in TelescopeFrame
Camera geometry, the cleaning mask should be applied to improve performance
image : array_like
Charge in each pixel, the cleaning mask should already be applied to
Expand Down Expand Up @@ -178,21 +178,6 @@ def hillas_parameters(geom, image):
np.sum(((((b * A) + (a * B) + (-c * C))) ** 2.0) * image)
) / (2 * width)

if unit.is_equivalent(u.m):
return CameraHillasParametersContainer(
x=u.Quantity(cog_x, unit),
y=u.Quantity(cog_y, unit),
r=u.Quantity(cog_r, unit),
phi=Angle(cog_phi, unit=u.rad),
intensity=size,
length=u.Quantity(length, unit),
length_uncertainty=u.Quantity(length_uncertainty, unit),
width=u.Quantity(width, unit),
width_uncertainty=u.Quantity(width_uncertainty, unit),
psi=Angle(psi, unit=u.rad),
skewness=skewness_long,
kurtosis=kurtosis_long,
)
return HillasParametersContainer(
fov_lon=u.Quantity(cog_x, unit),
fov_lat=u.Quantity(cog_y, unit),
Expand Down
47 changes: 12 additions & 35 deletions ctapipe/image/image_processor.py
Original file line number Diff line number Diff line change
@@ -1,23 +1,19 @@
"""
High level image processing (ImageProcessor Component)
"""
from copy import deepcopy

import numpy as np

from ctapipe.coordinates import TelescopeFrame

from ..containers import (
ArrayEventContainer,
CameraHillasParametersContainer,
CameraTimingParametersContainer,
ImageParametersContainer,
IntensityStatisticsContainer,
PeakTimeStatisticsContainer,
TimingParametersContainer,
)
from ..core import QualityQuery, TelescopeComponent
from ..core.traits import Bool, BoolTelescopeParameter, ComponentName, List
from ..core.traits import BoolTelescopeParameter, ComponentName, List
from ..instrument import SubarrayDescription
from .cleaning import ImageCleaner
from .concentration import concentration_parameters
Expand All @@ -40,15 +36,9 @@
kurtosis=np.float64(np.nan),
)
DEFAULT_TIMING_PARAMETERS = TimingParametersContainer()
DEFAULT_TIMING_PARAMETERS_CAMFRAME = CameraTimingParametersContainer()
DEFAULT_PEAKTIME_STATISTICS = PeakTimeStatisticsContainer()


DEFAULT_IMAGE_PARAMETERS_CAMFRAME = deepcopy(DEFAULT_IMAGE_PARAMETERS)
DEFAULT_IMAGE_PARAMETERS_CAMFRAME.hillas = CameraHillasParametersContainer()
DEFAULT_IMAGE_PARAMETERS_CAMFRAME.timing = CameraTimingParametersContainer()


class ImageQualityQuery(QualityQuery):
"""for configuring image-wise data checks"""

Expand All @@ -68,11 +58,6 @@ class ImageProcessor(TelescopeComponent):
ImageCleaner, default_value="TailcutsImageCleaner"
).tag(config=True)

use_telescope_frame = Bool(
default_value=True,
help="Whether to calculate parameters in the telescope or camera frame",
).tag(config=True)

apply_image_modifier = BoolTelescopeParameter(
default_value=False, help="If true, apply ImageModifier to dl1 images"
).tag(config=True)
Expand Down Expand Up @@ -105,16 +90,14 @@ def __init__(

self.check_image = ImageQualityQuery(parent=self)

self.default_image_container = DEFAULT_IMAGE_PARAMETERS_CAMFRAME
if self.use_telescope_frame:
self.default_image_container = DEFAULT_IMAGE_PARAMETERS
telescope_frame = TelescopeFrame()
self.telescope_frame_geometries = {
tel_id: self.subarray.tel[tel_id].camera.geometry.transform_to(
telescope_frame
)
for tel_id in self.subarray.tel
}
self.default_image_container = DEFAULT_IMAGE_PARAMETERS
telescope_frame = TelescopeFrame()
self.telescope_frame_geometries = {
tel_id: self.subarray.tel[tel_id].camera.geometry.transform_to(
telescope_frame
)
for tel_id in self.subarray.tel
}

def __call__(self, event: ArrayEventContainer):
self._process_telescope_event(event)
Expand Down Expand Up @@ -182,10 +165,7 @@ def _parameterize_image(
container_class=PeakTimeStatisticsContainer,
)
else:
if self.use_telescope_frame:
timing = DEFAULT_TIMING_PARAMETERS
else:
timing = DEFAULT_TIMING_PARAMETERS_CAMFRAME
timing = DEFAULT_TIMING_PARAMETERS

peak_time_statistics = DEFAULT_PEAKTIME_STATISTICS

Expand All @@ -209,11 +189,8 @@ def _process_telescope_event(self, event):
"""
for tel_id, dl1_camera in event.dl1.tel.items():

if self.use_telescope_frame:
# Use the transformed geometries
geometry = self.telescope_frame_geometries[tel_id]
else:
geometry = self.subarray.tel[tel_id].camera.geometry
# Use the transformed geometries
geometry = self.telescope_frame_geometries[tel_id]

if self.apply_image_modifier.tel[tel_id]:
dl1_camera.image = self.modify(tel_id=tel_id, image=dl1_camera.image)
Expand Down
11 changes: 6 additions & 5 deletions ctapipe/image/muon/tests/test_ring_fitter.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
import astropy.units as u
import pytest

from ctapipe.coordinates import TelescopeFrame
from ctapipe.image import tailcuts_clean, toymodel
from ctapipe.image.muon import MuonRingFitter

Expand All @@ -15,10 +16,10 @@ def test_MuonRingFitter_has_methods():
def test_MuonRingFitter(method, prod5_mst_flashcam):
"""test MuonRingFitter"""
# flashCam example
center_xs = 0.3 * u.m
center_ys = 0.6 * u.m
radius = 0.3 * u.m
width = 0.05 * u.m
center_xs = 0.3 * u.deg
center_ys = 0.6 * u.deg
radius = 0.5 * u.deg
width = 0.05 * u.deg

muon_model = toymodel.RingGaussian(
x=center_xs,
Expand All @@ -28,7 +29,7 @@ def test_MuonRingFitter(method, prod5_mst_flashcam):
)

# testing with flashcam
geom = prod5_mst_flashcam.camera.geometry
geom = prod5_mst_flashcam.camera.geometry.transform_to(TelescopeFrame())
charge, _, _ = muon_model.generate_image(
geom,
intensity=1000,
Expand Down
Loading