Skip to content

Commit

Permalink
Add EngineeringCameraFrame (#1085)
Browse files Browse the repository at this point in the history
* Add EngineeringCameraFrame

* Allow plotting in engineering frame

* Add transform_to to CameraGeometry and use it in plotting

* Rename constants

* Move transformation from plotting to Geometry

* Add example to the docs

* Remove unused imports

* Add dual-mirror case

* Style fixes

* Update docs

* Add missing docstring
  • Loading branch information
maxnoe authored Jun 17, 2019
1 parent c154718 commit 5eb9995
Show file tree
Hide file tree
Showing 8 changed files with 225 additions and 18 deletions.
3 changes: 2 additions & 1 deletion ctapipe/coordinates/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,13 @@
from .telescope_frame import TelescopeFrame
from .nominal_frame import NominalFrame
from .ground_frames import GroundFrame, TiltedGroundFrame, project_to_ground
from .camera_frame import CameraFrame
from .camera_frame import CameraFrame, EngineeringCameraFrame


__all__ = [
'TelescopeFrame',
'CameraFrame',
'EngineeringCameraFrame',
'NominalFrame',
'GroundFrame',
'TiltedGroundFrame',
Expand Down
73 changes: 73 additions & 0 deletions ctapipe/coordinates/camera_frame.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,19 +4,48 @@
BaseCoordinateFrame,
CoordinateAttribute,
QuantityAttribute,
Attribute,
TimeAttribute,
EarthLocationAttribute,
FunctionTransform,
frame_transform_graph,
CartesianRepresentation,
UnitSphericalRepresentation,
AltAz,
AffineTransform,
)

from .telescope_frame import TelescopeFrame
from .representation import PlanarRepresentation


class MirrorAttribute(Attribute):
'''A frame Attribute that can only store the integers 1 and 2'''

def convert_input(self, value):
'''make sure input is 1 or 2'''
if value in (1, 2):
return value, False

raise ValueError('Only 1 or 2 mirrors supported')


# Go from SimTel / HESS to MAGIC/FACT/Engineering frame and back
CAMERA_TO_ENGINEERING_1M_MATRIX = np.array([
[0, -1, 0],
[-1, 0, 0],
[0, 0, 1]
])
ENGINEERING_1M_TO_CAMERA_MATRIX = CAMERA_TO_ENGINEERING_1M_MATRIX
CAMERA_TO_ENGINEERING_2M_MATRIX = np.array([
[0, 1, 0],
[-1, 0, 0],
[0, 0, 1]
])
ENGINEERING_2M_TO_CAMERA_MATRIX = CAMERA_TO_ENGINEERING_2M_MATRIX.T
ZERO_OFFSET = CartesianRepresentation(0, 0, 0, unit=u.m)


class CameraFrame(BaseCoordinateFrame):
'''
Camera coordinate frame.
Expand Down Expand Up @@ -58,6 +87,35 @@ class CameraFrame(BaseCoordinateFrame):
location = EarthLocationAttribute(default=None)


class EngineeringCameraFrame(CameraFrame):
'''
Engineering camera coordinate frame.
The camera frame is a 2d cartesian frame,
describing position of objects in the focal plane of the telescope.
The frame is defined as in MAGIC and FACT.
Standing in the dish, looking onto the camera, x points right and y points up.
HESS, ctapipe and sim_telarray use a different camera coordinate system:
To transform H.E.S.S./ctapipe -> FACT/MAGIC, do x' = -y, y' = -x.
Attributes
----------
focal_length : u.Quantity[length]
Focal length of the telescope as a unit quantity (usually meters)
rotation : u.Quantity[angle]
Rotation angle of the camera (0 deg in most cases)
telescope_pointing : SkyCoord[AltAz]
Pointing direction of the telescope as SkyCoord in AltAz
obstime : Time
Observation time
location : EarthLocation
location of the telescope
'''
n_mirrors = MirrorAttribute(default=1)


@frame_transform_graph.transform(FunctionTransform, CameraFrame, TelescopeFrame)
def camera_to_telescope(camera_coord, telescope_frame):
'''
Expand Down Expand Up @@ -132,3 +190,18 @@ def telescope_to_camera(telescope_coord, camera_frame):
)

return camera_frame.realize_frame(representation)


@frame_transform_graph.transform(AffineTransform, CameraFrame, EngineeringCameraFrame)
def camera_to_engineering(from_coord, to_frame):
if to_frame.n_mirrors == 1:
return CAMERA_TO_ENGINEERING_1M_MATRIX, ZERO_OFFSET

return CAMERA_TO_ENGINEERING_2M_MATRIX, ZERO_OFFSET


@frame_transform_graph.transform(AffineTransform, EngineeringCameraFrame, CameraFrame)
def engineering_to_camera(from_coord, to_frame):
if from_coord.n_mirrors == 1:
return ENGINEERING_1M_TO_CAMERA_MATRIX, ZERO_OFFSET
return ENGINEERING_2M_TO_CAMERA_MATRIX, ZERO_OFFSET
33 changes: 33 additions & 0 deletions ctapipe/coordinates/tests/test_engineering_frame.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
'''
Tests for the conversion between camera coordinate frames
'''
from astropy.coordinates import SkyCoord
import astropy.units as u


def test_conversion():
'''
Test conversion between CameraFrame and EngineeringCameraFrame
'''
from ctapipe.coordinates import CameraFrame, EngineeringCameraFrame

coords = SkyCoord(x=[3, 1] * u.m, y=[2, 4] * u.m, frame=CameraFrame())

for coord in coords:
eng_coord = coord.transform_to(EngineeringCameraFrame())

assert eng_coord.x == -coord.y
assert eng_coord.y == -coord.x

back = eng_coord.transform_to(CameraFrame())
assert back.x == coord.x
assert back.y == coord.y

eng_coord = coord.transform_to(EngineeringCameraFrame(n_mirrors=2))

assert eng_coord.x == coord.y
assert eng_coord.y == -coord.x

back = eng_coord.transform_to(CameraFrame())
assert back.x == coord.x
assert back.y == coord.y
45 changes: 43 additions & 2 deletions ctapipe/instrument/camera.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@

import numpy as np
from astropy import units as u
from astropy.coordinates import Angle
from astropy.coordinates import Angle, SkyCoord
from astropy.table import Table
from astropy.utils import lazyproperty
from scipy.spatial import cKDTree as KDTree
Expand All @@ -15,6 +15,7 @@

from ctapipe.utils import get_table_dataset, find_all_matching_datasets
from ctapipe.utils.linalg import rotation_matrix_2d
from ctapipe.coordinates import CameraFrame


__all__ = ['CameraGeometry']
Expand Down Expand Up @@ -68,7 +69,7 @@ class CameraGeometry:

def __init__(self, cam_id, pix_id, pix_x, pix_y, pix_area, pix_type,
pix_rotation="0d", cam_rotation="0d",
neighbors=None, apply_derotation=True):
neighbors=None, apply_derotation=True, frame=None):

if pix_x.ndim != 1 or pix_y.ndim != 1:
raise ValueError(f'Pixel coordinates must be 1 dimensional, got {pix_x.ndim}')
Expand All @@ -84,6 +85,7 @@ def __init__(self, cam_id, pix_id, pix_x, pix_y, pix_area, pix_type,
self.pix_rotation = Angle(pix_rotation)
self.cam_rotation = Angle(cam_rotation)
self._neighbors = neighbors
self.frame = frame

if neighbors is not None:
if isinstance(neighbors, list):
Expand Down Expand Up @@ -125,6 +127,45 @@ def __eq__(self, other):
(self.pix_y == other.pix_y).all(),
])

def transform_to(self, frame):
'''
Transform the pixel coordinates stored in this geometry
and the pixel and camera rotations to another camera coordinate frame.
Parameters
----------
frame: ctapipe.coordinates.CameraFrame
The coordinate frame to transform to.
'''
if self.frame is None:
self.frame = CameraFrame()

coord = SkyCoord(x=self.pix_x, y=self.pix_y, frame=self.frame)
trans = coord.transform_to(frame)

# also transform the unit vectors, to get rotation / mirroring
uv = SkyCoord(x=[1, 0], y=[0, 1], unit=u.m, frame=self.frame)
uv_trans = uv.transform_to(frame)
rot = np.arctan2(uv_trans[0].y, uv_trans[1].y)
det = np.linalg.det([uv_trans.x.value, uv_trans.y.value])

cam_rotation = rot + det * self.cam_rotation
pix_rotation = rot + det * self.pix_rotation

return CameraGeometry(
cam_id=self.cam_id,
pix_id=self.pix_id,
pix_x=trans.x,
pix_y=trans.y,
pix_area=self.pix_area,
pix_type=self.pix_type,
pix_rotation=pix_rotation,
cam_rotation=cam_rotation,
neighbors=None,
apply_derotation=False,
frame=frame,
)

def __hash__(self):
return hash((
self.cam_id,
Expand Down
16 changes: 15 additions & 1 deletion ctapipe/instrument/tests/test_camera.py
Original file line number Diff line number Diff line change
Expand Up @@ -287,8 +287,22 @@ def test_hashing():

assert len(set([cam1, cam2, cam3])) == 2


@pytest.mark.parametrize("camera_name", CameraGeometry.get_known_camera_names())
def test_camera_from_name(camera_name):
""" check we can construct all cameras from name"""
camera = CameraGeometry.from_name(camera_name)
assert str(camera) == camera_name
assert str(camera) == camera_name


@pytest.mark.parametrize("camera_name", CameraGeometry.get_known_camera_names())
def test_camera_coordinate_transform(camera_name):
'''test conversion of the coordinates stored in a camera frame'''
from ctapipe.coordinates import EngineeringCameraFrame

geom = CameraGeometry.from_name(camera_name)
trans_geom = geom.transform_to(EngineeringCameraFrame())

unit = geom.pix_x.unit
assert np.allclose(geom.pix_x.to_value(unit), -trans_geom.pix_y.to_value(unit))
assert np.allclose(geom.pix_y.to_value(unit), -trans_geom.pix_x.to_value(unit))
29 changes: 15 additions & 14 deletions ctapipe/visualization/mpl_camera.py
Original file line number Diff line number Diff line change
Expand Up @@ -94,10 +94,9 @@ def __init__(
cmap=None,
allow_pick=False,
autoupdate=True,
autoscale=True
autoscale=True,
):
self.axes = ax if ax is not None else plt.gca()
self.geom = geometry
self.pixels = None
self.colorbar = None
self.autoupdate = autoupdate
Expand All @@ -106,6 +105,8 @@ def __init__(
self._active_pixel_label = None
self._axes_overlays = []

self.geom = geometry

if title is None:
title = geometry.cam_id

Expand All @@ -117,25 +118,25 @@ def __init__(
if not hasattr(self.geom, "mask"):
self.geom.mask = np.ones_like(self.geom.pix_x.value, dtype=bool)

for xx, yy, aa in zip(
u.Quantity(self.geom.pix_x[self.geom.mask]).value,
u.Quantity(self.geom.pix_y[self.geom.mask]).value,
u.Quantity(np.array(self.geom.pix_area)[self.geom.mask]).value):
pix_x = self.geom.pix_x.value[self.geom.mask]
pix_y = self.geom.pix_y.value[self.geom.mask]
pix_area = self.geom.pix_area.value[self.geom.mask]

for x, y, area in zip(pix_x, pix_y, pix_area):
if self.geom.pix_type.startswith("hex"):
rr = sqrt(aa * 2 / 3 / sqrt(3)) + 2 * PIXEL_EPSILON
r = sqrt(area * 2 / 3 / sqrt(3)) + 2 * PIXEL_EPSILON
poly = RegularPolygon(
(xx, yy), 6, radius=rr,
orientation=self.geom.pix_rotation.rad,
(x, y), 6, radius=r,
orientation=self.geom.pix_rotation.to_value(u.rad),
fill=True,
)
else:
rr = sqrt(aa) + PIXEL_EPSILON
r = sqrt(area) + PIXEL_EPSILON
poly = Rectangle(
(xx - rr / 2., yy - rr / 2.),
width=rr,
height=rr,
angle=self.geom.pix_rotation.deg,
(x - r / 2, y - r / 2),
width=r,
height=r,
angle=self.geom.pix_rotation.to_value(u.deg),
fill=True,
)

Expand Down
10 changes: 10 additions & 0 deletions docs/ctapipe_api/coordinates/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ Coordinates should be described using an `astropy.coordinates.SkyCoord` in the a
The following special frames are defined for CTA:

* `CameraFrame`
* `EngineeringCameraFrame`
* `TelescopeFrame`
* `NominalFrame`
* `GroundFrame`
Expand All @@ -35,6 +36,15 @@ The following special frames are defined for CTA:
they can be transformed to and from any other `astropy.coordinates` frame, like
`astropy.coordinates.AltAz` or `astropy.coordinates.ICRS` (RA/Dec)

The two different coordinate frames are shown here:

.. plot:: ../examples/plot_camera_frames.py
:include-source:

The `CameraFrame` is used internally in `ctapipe` and comes from `sim_telarray`.
It abstracts the transformation differences between 1 and 2 mirror telescopes away.
The `EngineeringCameraFrame` is used by `MAGIC`, `FACT` and the `H.E.S.S.` analysis
software.


Reference/API
Expand Down
34 changes: 34 additions & 0 deletions examples/plot_camera_frames.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
'''
Plot the same event in two camera displays showing the
different coordinate frames for camera coordinates.
'''
from ctapipe.instrument import CameraGeometry
from ctapipe.visualization import CameraDisplay
import matplotlib.pyplot as plt
from ctapipe.coordinates import EngineeringCameraFrame
from ctapipe.image.toymodel import Gaussian
import astropy.units as u


def main():
fig, axs = plt.subplots(1, 2, constrained_layout=True, figsize=(6, 3))

model = Gaussian(0 * u.m, 0.1 * u.m, 0.3 * u.m, 0.05 * u.m, 25 * u.deg)
cam = CameraGeometry.from_name('FlashCam')
image, *_ = model.generate_image(cam, 2500)

CameraDisplay(cam, ax=axs[0], image=image)
CameraDisplay(
cam.transform_to(EngineeringCameraFrame()),
ax=axs[1],
image=image,
)

axs[0].set_title('CameraFrame')
axs[1].set_title('EngineeringCameraFrame')

plt.show()


if __name__ == '__main__':
main()

0 comments on commit 5eb9995

Please sign in to comment.