diff --git a/README.rst b/README.rst
index 30ddda8421ec672dd282844d624ad899a1882028..34a9de28c89bcabb7db7e52e433897e06555a725 100644
--- a/README.rst
+++ b/README.rst
@@ -33,7 +33,8 @@ environment (this is optional).
The following command will set up a conda virtual environment, add the
necessary package channels, and download ctapipe and its dependencies. The
-file *environment.yml* can be found in this repo.
+file *environment.yml* can be found in this repo.
+Note this is *pre-alpha* software and is not yet stable enough for end-users (expect large API changes until the first stable 1.0 release).
::
diff --git a/ctapipe/coordinates/__init__.py b/ctapipe/coordinates/__init__.py
index 8bc3db6a1b171646efd096c1b3608dd028b1f85f..c5ba24ed797600d398296775804183c5bdb81382 100644
--- a/ctapipe/coordinates/__init__.py
+++ b/ctapipe/coordinates/__init__.py
@@ -4,9 +4,14 @@ Coordinates.
"""
-from astropy.utils import iers
-iers.conf.auto_download = True # auto-fetch updates to IERS_A table
+from .coordinate_transformations import (pixel_position_to_direction,
+ alt_to_theta, az_to_phi,
+ transform_pixel_position)
from .angular_frames import *
from .ground_frames import *
+
+
+from astropy.utils import iers
+iers.conf.auto_download = True # auto-fetch updates to IERS_A table
diff --git a/ctapipe/coordinates/coordinate_transformations.py b/ctapipe/coordinates/coordinate_transformations.py
new file mode 100644
index 0000000000000000000000000000000000000000..b5d3f0fb16f995d9d69735e02a114b696a8ec78e
--- /dev/null
+++ b/ctapipe/coordinates/coordinate_transformations.py
@@ -0,0 +1,82 @@
+import os
+
+import numpy as np
+from astropy.utils.decorators import deprecated
+from astropy import units as u
+
+from ctapipe import utils
+from ctapipe.utils import linalg
+
+
+@deprecated(0.1, "will be replaced with proper coord transform")
+def pixel_position_to_direction(pix_x, pix_y, tel_phi, tel_theta, tel_foclen):
+ """
+ TODO replace with proper implementation
+ calculates the direction vector of corresponding to a
+ (x,y) position on the camera
+
+ beta is the pixel's angular distance to the centre
+ according to beta / tel_view = r / maxR
+ alpha is the polar angle between the y-axis and the pixel
+ to find the direction the pixel is looking at:
+
+ - the pixel direction is set to the telescope direction
+ - offset by beta towards up
+ - rotated around the telescope direction by the angle alpha
+
+
+ Parameters
+ -----------
+ pix_x, pix_y : ndarray
+ lists of x and y positions on the camera
+ tel_phi, tel_theta: astropy quantities
+ two angles that describe the orientation of the telescope
+ tel_foclen : astropy quantity
+ focal length of the telescope
+
+ Returns
+ -------
+ pix_dirs : ndarray
+ shape (n,3) list of "direction vectors"
+ corresponding to a position on the camera
+ """
+
+ pix_alpha = np.arctan2(pix_y, pix_x)
+
+ pix_rho = (pix_x ** 2 + pix_y ** 2) ** .5
+
+ pix_beta = pix_rho / tel_foclen * u.rad
+
+ tel_dir = linalg.set_phi_theta(tel_phi, tel_theta)
+
+ pix_dirs = []
+ for a, b in zip(pix_alpha, pix_beta):
+ pix_dir = linalg.set_phi_theta(tel_phi, tel_theta - b)
+
+ pix_dir = linalg.rotate_around_axis(pix_dir, tel_dir, a)
+ pix_dirs.append(pix_dir * u.dimless)
+
+ return pix_dirs
+
+
+def alt_to_theta(alt):
+ """transforms altitude (angle from the horizon upwards) to theta (angle from z-axis)
+ for simtel array coordinate systems
+ """
+ return (90 * u.deg - alt).to(alt.unit)
+
+
+def az_to_phi(az):
+ """transforms azimuth (angle from north towards east)
+ to phi (angle from x-axis towards y-axis)
+ for simtel array coordinate systems
+ """
+ return -az
+
+
+def transform_pixel_position(x, y):
+ """transforms the x and y coordinates on the camera plane so that they correspond to
+ the view as if looked along the pointing direction of the telescope, i.e. y->up and
+ x->right
+ """
+ return x, -y
diff --git a/ctapipe/core/component.py b/ctapipe/core/component.py
index 6badade143e9ded7d1fa5e2c1ce3bffe07c12fdb..ddfa3d6158e4c0751957de6d706904a914463561 100644
--- a/ctapipe/core/component.py
+++ b/ctapipe/core/component.py
@@ -15,7 +15,7 @@ class AbstractConfigurableMeta(type(Configurable), ABCMeta):
class Component(Configurable, metaclass=AbstractConfigurableMeta):
"""Base class of all Components (sometimes called
- workers, makers, etc). Components are are classes that do some sort
+ workers, makers, etc). Components are classes that do some sort
of processing and contain user-configurable parameters, which are
implemented using `traitlets`.
@@ -61,7 +61,7 @@ class Component(Configurable, metaclass=AbstractConfigurableMeta):
parent: Tool or Component
Tool or component that is the Parent of this one
kwargs: type
- other paremeters
+ other parameters
"""
@@ -71,4 +71,6 @@ class Component(Configurable, metaclass=AbstractConfigurableMeta):
if self.parent:
self.log = self.parent.log.getChild(self.__class__.__name__)
else:
- self.log = getLogger(self.__class__.__name__)
+ self.log = getLogger(
+ self.__class__.__module__ + '.' + self.__class__.__name__
+ )
diff --git a/ctapipe/core/provenance.py b/ctapipe/core/provenance.py
index 83a805e2fde3763bec4d54e756477516c069d6ab..fbfd5a6626ea0012d53d7c9cc3c979209ff66736 100644
--- a/ctapipe/core/provenance.py
+++ b/ctapipe/core/provenance.py
@@ -147,7 +147,7 @@ class Provenance(metaclass=Singleton):
class _ActivityProvenance:
"""
Low-level helper class to collect provenance information for a given
- *activity*. Users should use `Provenance` as a top-level API,
+ *activity*. Users should use `Provenance` as a top-level API,
not this class directly.
"""
@@ -180,9 +180,9 @@ class _ActivityProvenance:
url: str
filename or url of input file
role: str
- role name that this output satisfies
+ role name that this input satisfies
"""
- self._prov['input'].append(dict(url=url,role=role))
+ self._prov['input'].append(dict(url=url, role=role))
def register_output(self, url, role=None):
"""
@@ -196,7 +196,7 @@ class _ActivityProvenance:
role: str
role name that this output satisfies
"""
- self._prov['output'].append(dict(url=url,role=role))
+ self._prov['output'].append(dict(url=url, role=role))
def register_config(self, config):
""" add a dictionary of configuration parameters to this activity"""
diff --git a/ctapipe/core/tool.py b/ctapipe/core/tool.py
index 10b18a077c3c639b78b82e936cacbd02b2b80c1a..71c85573293e13bfe47d396479ef8bb2fd25cb9e 100644
--- a/ctapipe/core/tool.py
+++ b/ctapipe/core/tool.py
@@ -4,12 +4,12 @@ from abc import abstractmethod
from traitlets import Unicode
from traitlets.config import Application
-logging.basicConfig(level=logging.WARNING)
-
from ctapipe import __version__ as version
from .logging import ColoredFormatter
from . import Provenance
+logging.basicConfig(level=logging.WARNING)
+
class ToolConfigurationError(Exception):
@@ -113,9 +113,9 @@ class Tool(Application):
self.aliases['config'] = 'Tool.config_file'
super().__init__(**kwargs)
- self.log_format = ('%(levelname)8s [%(name)s] '
+ self.log_format = ('%(levelname)8s [%(name)s] '
'(%(module)s/%(funcName)s): %(message)s')
- self.log_level = 20 # default to INFO and above
+ self.log_level = logging.INFO
self.is_setup = False
def initialize(self, argv=None):
@@ -153,7 +153,7 @@ class Tool(Application):
----------
argv: list(str)
- command-line arguments, or None to get them
+ command-line arguments, or None to get them
from sys.argv automatically
"""
try:
@@ -182,7 +182,7 @@ class Tool(Application):
status='interrupted')
finally:
for activity in Provenance().finished_activities:
- output_str = ' '.join([x['url'] for x in activity.output])
+ output_str = ' '.join([x['url'] for x in activity.output])
self.log.info("Output: %s", output_str)
self.log.debug("PROVENANCE: '%s'", Provenance().as_json(indent=3))
diff --git a/ctapipe/flow/flow.py b/ctapipe/flow/flow.py
index 85f6e99d91fb82736d56a00f29ad7875ee6bf7c8..d88696a56baacff38d5865ceec51482921e0f1ec 100644
--- a/ctapipe/flow/flow.py
+++ b/ctapipe/flow/flow.py
@@ -31,7 +31,7 @@ class PipeStep():
'''
PipeStep represents a Flow step. One or several processes can be attach
- to this step.
+ to this step.
Parameters
----------
name : str
@@ -73,14 +73,15 @@ class PipeStep():
def __repr__(self):
'''standard representation
'''
- return ('Name[ ' + str(self.name)
- + '], next_steps_name[' + str(self.next_steps_name)
- + '], port in[ ' + str(self.port_in)
- + '], main connection name [ ' + str(self.main_connection_name) + ' ]'
- + '], port in[ ' + str(self.port_in)
- + '], nb process[ ' + str(self.nb_process)
- + '], level[ ' + str(self.level)
- + '], queue_limit[ ' + str(self.queue_limit) + ']')
+ return (
+ 'Name[{self.name}], '
+ 'next_steps_name[{self.next_steps_name}], '
+ 'port in[{self.port_in}], '
+ 'main connection name [{self.main_connection_name} ], '
+ 'nb_process [{self.nb_process}], '
+ 'level [{self.level}]'
+ 'queue_limit [{self.queue_limit}]'
+ ).format(self=self)
class FlowError(Exception):
@@ -269,11 +270,11 @@ class Flow(Tool):
return True
def configure_stagers(self, router_names):
- """ Creates Processes with users's coroutines for all stages
+ """ Creates Processes with users's coroutines for all stages
Parameters
----------
router_names: List
- List to fill with routers name
+ List to fill with routers name
Returns
-------
True if every instantiation is correct
@@ -310,7 +311,7 @@ class Flow(Tool):
def configure_consumer(self):
- """ Creates consumer Processes with users's coroutines
+ """ Creates consumer Processes with users's coroutines
Returns
-------
True if every instantiation is correct
@@ -329,7 +330,7 @@ class Flow(Tool):
def add_consumer_to_router(self):
""" Create router_names dictionary and
- Add consumer router ports
+ Add consumer router ports
Returns
-------
The new router_names dictionary
@@ -344,7 +345,7 @@ class Flow(Tool):
return router_names
def configure_producer(self):
- """ Creates producer Process with users's coroutines
+ """ Creates producer Process with users's coroutines
Returns
-------
True if every instatiation is correct
@@ -365,7 +366,7 @@ class Flow(Tool):
return True
def connect_gui(self):
- """ Connect ZMQ socket to send information to GUI
+ """ Connect ZMQ socket to send information to GUI
Returns
-------
True if everything correct
@@ -402,7 +403,7 @@ class Flow(Tool):
def configure_ports(self):
"""
- Configures producer, stagers and consumer ZMQ ports
+ Configures producer, stagers and consumer ZMQ ports
Returns
-------
True if everything correct
@@ -444,11 +445,11 @@ class Flow(Tool):
def get_step_by_name(self, name):
''' Find a PipeStep in self.producer_step or self.stager_steps or
- self.consumer_step
+ self.consumer_step
Parameters
----------
name : str
- step name
+ step name
Returns
-------
PipeStep if found, otherwise None
@@ -461,7 +462,7 @@ class Flow(Tool):
def instantiation(self, name, stage_type, process_name=None, port_in=None,
connections=None, main_connection_name=None, config=None):
'''
- Instantiate on Python object from name found in configuration
+ Instantiate on Python object from name found in configuration
Parameters
----------
name : str
@@ -504,13 +505,13 @@ class Flow(Tool):
def get_pipe_steps(self, role):
'''
- Create a list of Flow based framework steps from configuration and
- filter by role
+ Create a list of Flow based framework steps from configuration and
+ filter by role
Parameters
----------
role: str
filter with role for step to be add in result list
- Accepted values: self.PRODUCER - self.STAGER - self.CONSUMER
+ Accepted values: self.PRODUCER - self.STAGER - self.CONSUMER
Returns
-------
PRODUCER,CONSUMER: a step name filter by specific role (PRODUCER,CONSUMER)
@@ -557,7 +558,7 @@ class Flow(Tool):
return None
def def_step_for_gui(self):
- '''
+ '''
Create a list (levels_for_gui) containing all steps
Returns
@@ -586,7 +587,7 @@ class Flow(Tool):
for process in step.process:
nb_job_done += process.nb_job_done
running += process.running
- levels_for_gui.append(StagerRep(process.name, step.next_steps_name,
+ levels_for_gui.append(StagerRep(process.name, step.next_steps_name,
nb_job_done=nb_job_done,
running=running,
nb_process=len(step.process)))
@@ -680,12 +681,12 @@ class Flow(Tool):
def run_generator(self, destination, msg):
""" Get step for destination. Create a genetor from its run method.
- re-enter in run_generator until Generator send values
+ re-enter in run_generator until Generator send values
Parameters
----------
destination: str
Next step name
- msg: a Pickle dumped msg
+ msg: a Pickle dumped msg
Returns
-------
Next destination and msg
@@ -769,7 +770,7 @@ class Flow(Tool):
def wait_all_stagers(self, mintime):
""" Verify id all steps (stage + consumers) are finished their
- jobs and waiting
+ jobs and waiting
Returns
-------
True if all stages queue are empty and all Processes
@@ -789,9 +790,9 @@ class Flow(Tool):
def wait_and_send_levels(self, processes_to_wait):
'''
- Wait for a process to join and regularly send Flow based framework
+ Wait for a process to join and regularly send Flow based framework
state to GUI
- in case of a GUI will connect later
+ in case of a GUI will connect later
Parameters
----------
processes_to_wait : process
@@ -811,7 +812,7 @@ class Flow(Tool):
def get_step_conf(self, name):
'''
Search step by its name in self.stage_conf list,
- self.producer_conf and self.consumer_conf
+ self.producer_conf and self.consumer_conf
Parameters
----------
name : str
@@ -832,11 +833,11 @@ class Flow(Tool):
def get_stager_indice(self, name):
'''
- Search step by its name in self.stage_conf list
+ Search step by its name in self.stage_conf list
Parameters
----------
name : str
- stage name
+ stage name
Returns
-------
indice in list, -1 if not found
diff --git a/ctapipe/image/geometry_converter.py b/ctapipe/image/geometry_converter.py
index 885e9993363ac5b9fc2898c6d423fe2785a4abff..e9b8b4adae4c574f28eec2fd8fafbd15e5100da8 100644
--- a/ctapipe/image/geometry_converter.py
+++ b/ctapipe/image/geometry_converter.py
@@ -1,705 +1,9 @@
-import logging
-from collections import namedtuple
-import numpy as np
-from astropy import units as u
-from numba import jit
+"""collects converter functions from the various `geometry_converter_XXX` files in one
+common module"""
-from ctapipe.instrument import CameraGeometry
-logger = logging.getLogger(__name__)
+from .geometry_converter_hex import (convert_geometry_hex1d_to_rect2d,
+ convert_geometry_rect2d_back_to_hexe1d)
-RotBuffer = namedtuple("RotBuffer",
- "rot_x,rot_y,x_edges,y_edges,new_geom,rot_angle,pix_rotation,"
- "x_scale")
-
-
-__all__ = [
- "convert_geometry_hex1d_to_rect2d",
- "convert_geometry_rect2d_back_to_hexe1d"
-]
-
-
-def unskew_hex_pixel_grid(pix_x, pix_y, cam_angle=0 * u.deg,
- base_angle=60 * u.deg):
- r"""transform the pixel coordinates of a hexagonal image into an
- orthogonal image
-
- Parameters
- ----------
- pix_x, pix_y : 1D numpy arrays
- the list of x and y coordinates of the hexagonal pixel grid
- cam_angle : astropy.Quantity (default: 0 degrees)
- The skewing is performed along the y-axis, therefore, one of the slanted
- base-vectors needs to be parallel to the y-axis.
- Some camera grids are rotated in a way that this is not the case.
- This needs to be corrected.
- base_angle : astropy.Quantity (default: 60 degrees)
- the skewing angle of the hex-grid. should be 60° for regular hexagons
-
- Returns
- -------
- pix_x, pix_y : 1D numpy arrays
- the list of x and y coordinates of the slanted, orthogonal pixel grid
-
- Notes
- -----
- The correction on the pixel position r can be described by a rotation R around
- one angle and a sheer S along a certain axis:
-
- .. math::
- r' = S \cdot R \cdot r
-
- .. math::
- \begin{pmatrix}
- x' \\
- y'
- \end{pmatrix}
- =
- \begin{pmatrix}
- 1 & 0 \\
- -1/\tan & 1
- \end{pmatrix}
- \cdot
- \begin{pmatrix}
- \cos & -\sin \\
- \sin & \cos
- \end{pmatrix}
- \cdot
- \begin{pmatrix}
- x \\
- y
- \end{pmatrix}
-
- .. math::
- \begin{pmatrix}
- x' \\
- y'
- \end{pmatrix}
- =
- \begin{pmatrix}
- \cos & -\sin \\
- \sin-\cos/\tan & \sin/\tan+\cos
- \end{pmatrix}
- \cdot
- \begin{pmatrix}
- x \\
- y
- \end{pmatrix}
-
- """
-
- tan_angle = np.tan(base_angle)
-
- # If the camera-rotation angle is non-zero, create a rotation+sheering
- # matrix for the pixel coordinates
- if cam_angle != 0 * u.deg:
- sin_angle = np.sin(cam_angle)
- cos_angle = np.cos(cam_angle)
-
- # the correction on the pixel position r can be described by a
- # rotation R around one angle and a sheer S along a certain axis:
- #
- # r' = S * R * r
- # (x') = ( 1 0) * (cos -sin) * (x) = ( cos -sin ) * (x)
- # (y') (-1/tan 1) (sin cos) (y) (sin-cos/tan sin/tan+cos) * (y)
- rot_mat = np.array(
- [[cos_angle, -sin_angle],
- [sin_angle - cos_angle / tan_angle,
- sin_angle / tan_angle + cos_angle]])
-
- else:
- # if we don't rotate the camera, only perform the sheer
- rot_mat = np.array([[1, 0], [-1 / tan_angle, 1]])
-
- rotated = np.dot(rot_mat, [pix_x.value, pix_y.value])
- rot_x = rotated[0] * pix_x.unit
- rot_y = rotated[1] * pix_x.unit
- return rot_x, rot_y
-
-
-def reskew_hex_pixel_grid(pix_x, pix_y, cam_angle=0 * u.deg,
- base_angle=60 * u.deg):
- r"""skews the orthogonal coordinates back to the hexagonal ones
-
- Parameters
- ----------
- pix_x, pix_y : 1D numpy arrays
- the list of x and y coordinates of the slanted, orthogonal pixel grid
- cam_angle : astropy.Quantity (default: 0 degrees)
- The skewing is performed along the y-axis, therefore, one of the slanted
- base-vectors needs to be parallel to the y-axis.
- Some camera grids are rotated in a way that this is not the case.
- This needs to be corrected.
- base_angle : astropy.Quantity (default: 60 degrees)
- the skewing angle of the hex-grid. should be 60° for regular hexagons
-
- Returns
- -------
- pix_x, pix_y : 1D numpy arrays
- the list of x and y coordinates of the hexagonal pixel grid
-
- Notes
- -----
- To revert the rotation, we need to find matrices S' and R' with
- :math:`S' \cdot S = 1` and :math:`R' \cdot R = 1`,
- so that :math:`r = R' \cdot S' \cdot S \cdot R \cdot r = R' \cdot S' \cdot r'`:
-
- .. math::
- \begin{pmatrix}
- x \\
- y
- \end{pmatrix}
- =
- \begin{pmatrix}
- \cos & \sin \\
- -\sin & \cos
- \end{pmatrix}
- \cdot
- \begin{pmatrix}
- 1 & 0 \\
- 1/\tan & 1
- \end{pmatrix}
- \cdot
- \begin{pmatrix}
- x' \\
- y'
- \end{pmatrix}
-
- .. math::
- \begin{pmatrix}
- x \\
- y
- \end{pmatrix}
- =
- \begin{pmatrix}
- \cos+\sin/\tan & \sin \\
- \cos/\tan-\sin & \cos
- \end{pmatrix}
- \cdot
- \begin{pmatrix}
- x' \\
- y'
- \end{pmatrix}
-
- """
-
- tan_angle = np.tan(base_angle)
-
- # If the camera-rotation angle is non-zero, create a rotation+sheering
- # matrix for the pixel coordinates
- if cam_angle != 0 * u.deg:
- sin_angle = np.sin(cam_angle)
- cos_angle = np.cos(cam_angle)
-
- # to revert the rotation, we need to find matrices S' and R'
- # S' * S = 1 and R' * R = 1
- # so that
- # r = R' * S' * S * R * r = R' * S' * r'
- #
- # (x) = ( cos sin) * ( 1 0) * (x') = (cos+sin/tan sin) * (x')
- # (y) (-sin cos) (1/tan 1) (y') (cos/tan-sin cos) (y')
-
- rot_mat = np.array(
- [[cos_angle + sin_angle / tan_angle, sin_angle],
- [cos_angle / tan_angle - sin_angle, cos_angle]])
-
- else:
- # if we don't rotate the camera, only perform the sheer
- rot_mat = np.array([[1, 0], [1 / tan_angle, 1]])
-
- rotated = np.dot(rot_mat, [pix_x.value, pix_y.value])
- rot_x = rotated[0] * pix_x.unit
- rot_y = rotated[1] * pix_x.unit
- return rot_x, rot_y
-
-
-@jit
-def reskew_hex_pixel_from_orthogonal_edges(x_edges, y_edges, square_mask):
- """extracts and skews the pixel coordinates from a 2D orthogonal
- histogram (i.e. the bin-edges) and skews them into the hexagonal
- image while selecting only the pixel that are selected by the
- given mask
-
- Parameters
- ----------
- x_edges, y_edges : 1darrays
- the bin edges of the 2D histogram
- square_mask : 2darray
- mask that selects the pixels actually belonging to the camera
-
- Returns
- -------
- unrot_x, unrot_y : 1darrays
- pixel coordinated reskewed into the hexagonal camera grid
- """
-
- unrot_x, unrot_y = [], []
- for i, x in enumerate((x_edges[:-1] + x_edges[1:]) / 2):
- for j, y in enumerate((y_edges[:-1] + y_edges[1:]) / 2):
- if square_mask[i][j]:
- x_unrot, y_unrot = reskew_hex_pixel_grid(x, y)
- unrot_x.append(x_unrot)
- unrot_y.append(y_unrot)
- return unrot_x, unrot_y
-
-
-@jit
-def get_orthogonal_grid_edges(pix_x, pix_y, scale_aspect=True):
- """calculate the bin edges of the slanted, orthogonal pixel grid to
- resample the pixel signals with np.histogramdd right after.
-
- Parameters
- ----------
- pix_x, pix_y : 1D numpy arrays
- the list of x and y coordinates of the slanted, orthogonal pixel grid
- scale_aspect : boolean (default: True)
- if True, rescales the x-coordinates to create square pixels
- (instead of rectangular ones)
-
- Returns
- --------
- x_edges, y_edges : 1D numpy arrays
- the bin edges for the slanted, orthogonal pixel grid
- x_scale : float
- factor by which the x-coordinates have been scaled
- """
-
- # finding the size of the square patches
- d_x = 99 * u.m
- d_y = 99 * u.m
- x_base = pix_x[0]
- y_base = pix_y[0]
- for x, y in zip(pix_x, pix_y):
- if abs(y - y_base) < abs(x - x_base):
- d_x = min(d_x, abs(x - x_base))
- for x, y in zip(pix_x, pix_y):
- if abs(y - y_base) > abs(x - x_base):
- d_y = min(d_y, abs(y - y_base))
-
- x_scale = 1
- if scale_aspect:
- x_scale = d_y / d_x
- pix_x *= x_scale
- d_x = d_y
-
- # with the maximal extension of the axes and the size of the pixels,
- # determine the number of bins in each direction
- NBinsx = np.around(abs(max(pix_x) - min(pix_x)) / d_x) + 2
- NBinsy = np.around(abs(max(pix_y) - min(pix_y)) / d_y) + 2
- x_edges = np.linspace(min(pix_x).value, max(pix_x).value, NBinsx)
- y_edges = np.linspace(min(pix_y).value, max(pix_y).value, NBinsy)
-
- return x_edges, y_edges, x_scale
-
-
-# add_angle = 180 * u.deg
-rot_buffer = {}
-
-
-def convert_geometry_1d_to_2d(geom, signal, key=None, add_rot=0):
- """converts the geometry object of a camera with a hexagonal grid into
- a square grid by slanting and stretching the 1D arrays of pixel x
- and y positions and signal intensities are converted to 2D
- arrays. If the signal array contains a time-dimension it is
- conserved.
-
- Parameters
- ----------
- geom : CameraGeometry object
- geometry object of hexagonal cameras
- signal : ndarray
- 1D (no timing) or 2D (with timing) array of the pmt signals
- key : (default: None)
- arbitrary key to store the transformed geometry in a buffer
- add_rot : int/float (default: 0)
- parameter to apply an additional rotation of @add_rot times 60°
-
- Returns
- -------
- new_geom : CameraGeometry object
- geometry object of the slanted picture now with a rectangular
- grid and a 2D grid for the pixel positions contains now a 2D
- masking array signifying which of the pixels came from the
- original geometry and which are simply fillers from the
- rectangular grid square_img : ndarray 2D (no timing) or 3D
- (with timing) array of the pmt signals
-
- """
-
- if key in rot_buffer:
-
- # if the conversion with this key was done and stored before,
- # just read it in
- (rot_x, rot_y, x_edges, y_edges, new_geom,
- rot_angle, pix_rotation, x_scale) = rot_buffer[key]
- else:
-
- # otherwise, we have to do the conversion now first, skew all
- # the coordinates of the original geometry
-
- # extra_rot is the angle to get back to aligned hexagons with flat
- # tops. Note that the pixel rotation angle brings the camera so that
- # hexagons have a point at the top, so need to go 30deg back to
- # make them flat
- extra_rot = geom.pix_rotation - 30 * u.deg
-
- # total rotation angle:
- rot_angle = (add_rot * 60 * u.deg) - extra_rot
- # if geom.cam_id.startswith("NectarCam")\
- # or geom.cam_id.startswith("LSTCam"):
- # rot_angle += geom.cam_rotation + 90 * u.deg
-
- logger.debug("geom={}".format(geom))
- logger.debug("rot={}, extra={}".format(rot_angle, extra_rot))
-
- rot_x, rot_y = unskew_hex_pixel_grid(geom.pix_x, geom.pix_y,
- cam_angle=rot_angle)
-
- # with all the coordinate points, we can define the bin edges
- # of a 2D histogram
- x_edges, y_edges, x_scale = get_orthogonal_grid_edges(rot_x, rot_y)
-
- # this histogram will introduce bins that do not correspond to
- # any pixel from the original geometry. so we create a mask to
- # remember the true camera pixels by simply throwing all pixel
- # positions into numpy.histogramdd: proper pixels contain the
- # value 1, false pixels the value 0.
- square_mask = np.histogramdd([rot_y, rot_x],
- bins=(y_edges, x_edges))[0].astype(bool)
-
- # to be consistent with the pixel intensity, instead of saving
- # only the rotated positions of the true pixels (rot_x and
- # rot_y), create 2D arrays of all x and y positions (also the
- # false ones).
- grid_x, grid_y = np.meshgrid((x_edges[:-1] + x_edges[1:]) / 2.,
- (y_edges[:-1] + y_edges[1:]) / 2.)
-
- ids = []
- # instead of blindly enumerating all pixels, let's instead
- # store a list of all valid -- i.e. picked by the mask -- 2D
- # indices
- for i, row in enumerate(square_mask):
- for j, val in enumerate(row):
- if val is True:
- ids.append((i, j))
-
- # the area of the pixels (note that this is still a deformed
- # image)
- pix_area = np.ones_like(grid_x) \
- * (x_edges[1] - x_edges[0]) * (y_edges[1] - y_edges[0])
-
- # creating a new geometry object with the attributes we just determined
- new_geom = CameraGeometry(
- cam_id=geom.cam_id + "_rect",
- pix_id=ids, # this is a list of all the valid coordinate pairs now
- pix_x=grid_x * u.m,
- pix_y=grid_y * u.m,
- pix_area=pix_area * u.m ** 2,
- neighbors=geom.neighbors,
- pix_type='rectangular', apply_derotation=False)
-
- # storing the pixel mask and camera rotation for later use
- new_geom.mask = square_mask
-
- if key is not None:
- # if a key is given, store the essential objects in a buffer
- rot_buffer[key] = RotBuffer(rot_x, rot_y, x_edges, y_edges, new_geom,
- rot_angle, geom.pix_rotation, x_scale)
-
- # resample the signal array to correspond to the square grid --
- # for signal arrays containing time slices (ndim > 1) or not
- # approach is the same as used for the mask only with the signal
- # as bin-weights
- if signal.ndim > 1:
- t_dim = signal.shape[1]
- square_img = np.histogramdd([np.repeat(rot_y, t_dim),
- np.repeat(rot_x, t_dim),
- [a for a in range(t_dim)] * len(rot_x)],
- bins=(y_edges, x_edges, range(t_dim + 1)),
- weights=signal.ravel())[0]
- else:
- square_img = np.histogramdd([rot_y, rot_x],
- bins=(y_edges, x_edges),
- weights=signal)[0]
-
- return new_geom, square_img
-
-
-unrot_buffer = {} # todo: should just make this a function with @lru_cache
-
-
-def convert_geometry_back(geom, signal, key, add_rot=0):
- """reverts the geometry distortion performed by
- convert_geometry_1d_to_2d back to a hexagonal grid stored in 1D
- arrays
-
- Parameters
- ----------
- geom : CameraGeometry
- geometry object where pixel positions are stored in a 2D
- rectangular camera grid
- signal : ndarray
- pixel intensity stored in a 2D rectangular camera grid
- key:
- key to retrieve buffered geometry information
- add_rot : int/float (default: 0)
- parameter to apply an additional rotation of @add_rot times 60°
-
- Returns
- -------
- unrot_geom : CameraGeometry
- pixel rotated back to a hexagonal grid stored in a 1D array
- signal : ndarray
- 1D (no timing) or 2D (with timing) array of the pmt signals
-
- """
- global unrot_buffer
-
- square_mask = geom.mask
-
- if key in unrot_buffer:
- unrot_geom = unrot_buffer[key]
- else:
- if key in rot_buffer:
- x_scale = rot_buffer[key].x_scale
- rot_angle = rot_buffer[key].rot_angle
- pix_rotation = rot_buffer[key].pix_rotation
- else:
- raise KeyError("key '{}' not found in the buffer".format(key))
-
- grid_x, grid_y = geom.pix_x / x_scale, geom.pix_y
-
- unrot_x, unrot_y = reskew_hex_pixel_grid(grid_x[square_mask],
- grid_y[square_mask],
- rot_angle)
-
- # TODO: probably should use base constructor, not guess here:
- # unrot_geom = CameraGeometry.guess(unrot_x, unrot_y, foc_len,
- # apply_derotation=False)
- unrot_geom = CameraGeometry(cam_id=geom.cam_id + "_hex",
- pix_id=np.arange(len(unrot_x)),
- pix_x=unrot_x,
- pix_y=unrot_y,
- pix_area=None, # recalc
- pix_type='hexagonal',
- pix_rotation=pix_rotation,
- # apply_derotation=False
- )
-
- unrot_buffer[key] = unrot_geom
-
- return unrot_geom, signal[square_mask, ...]
-
-
-def convert_geometry_hex1d_to_rect2d(geom, signal, key=None, add_rot=0):
- """converts the geometry object of a camera with a hexagonal grid into
- a square grid by slanting and stretching the 1D arrays of pixel x
- and y positions and signal intensities are converted to 2D
- arrays. If the signal array contains a time-dimension it is
- conserved.
-
- Parameters
- ----------
- geom : CameraGeometry object
- geometry object of hexagonal cameras
- signal : ndarray
- 1D (no timing) or 2D (with timing) array of the pmt signals
- key : (default: None)
- arbitrary key to store the transformed geometry in a buffer
- add_rot : int/float (default: 0)
- parameter to apply an additional rotation of `add_rot` times 60°
-
- Returns
- -------
- new_geom : CameraGeometry object
- geometry object of the slanted picture now with a rectangular
- grid and a 2D grid for the pixel positions. contains now a 2D
- masking array signifying which of the pixels came from the
- original geometry and which are simply fillers from the
- rectangular grid
- rot_img : ndarray 2D (no timing) or 3D (with timing)
- the rectangular signal image
- """
-
- if key in rot_buffer:
-
- # if the conversion with this key was done before and stored,
- # just read it in
- (geom, new_geom, hex_to_rect_map) = rot_buffer[key]
- else:
-
- # otherwise, we have to do the conversion first now,
- # skew all the coordinates of the original geometry
-
- # extra_rot is the angle to get back to aligned hexagons with flat
- # tops. Note that the pixel rotation angle brings the camera so that
- # hexagons have a point at the top, so need to go 30deg back to
- # make them flat
- extra_rot = geom.pix_rotation - 30 * u.deg
-
- # total rotation angle:
- rot_angle = (add_rot * 60 * u.deg) - extra_rot
-
- logger.debug("geom={}".format(geom))
- logger.debug("rot={}, extra={}".format(rot_angle, extra_rot))
-
- rot_x, rot_y = unskew_hex_pixel_grid(geom.pix_x, geom.pix_y,
- cam_angle=rot_angle)
-
- # with all the coordinate points, we can define the bin edges
- # of a 2D histogram
- x_edges, y_edges, x_scale = get_orthogonal_grid_edges(rot_x, rot_y)
-
- # this histogram will introduce bins that do not correspond to
- # any pixel from the original geometry. so we create a mask to
- # remember the true camera pixels by simply throwing all pixel
- # positions into numpy.histogramdd: proper pixels contain the
- # value 1, false pixels the value 0.
- square_mask = np.histogramdd([rot_y, rot_x],
- bins=(y_edges, x_edges))[0].astype(bool)
-
- # to be consistent with the pixel intensity, instead of saving
- # only the rotated positions of the true pixels (rot_x and
- # rot_y), create 2D arrays of all x and y positions (also the
- # false ones).
- grid_x, grid_y = np.meshgrid((x_edges[:-1] + x_edges[1:]) / 2.,
- (y_edges[:-1] + y_edges[1:]) / 2.)
-
- ids = []
- # instead of blindly enumerating all pixels, let's instead
- # store a list of all valid -- i.e. picked by the mask -- 2D
- # indices
- for i, row in enumerate(square_mask):
- for j, val in enumerate(row):
- if val is True:
- ids.append((i, j))
-
- # the area of the pixels (note that this is still a deformed
- # image)
- pix_area = np.ones_like(grid_x) \
- * (x_edges[1] - x_edges[0]) * (y_edges[1] - y_edges[0])
-
- # creating a new geometry object with the attributes we just determined
- new_geom = CameraGeometry(
- cam_id=geom.cam_id + "_rect",
- pix_id=ids, # this is a list of all the valid coordinate pairs now
- pix_x=grid_x * u.m,
- pix_y=grid_y * u.m,
- pix_area=pix_area * u.m ** 2,
- neighbors=geom.neighbors,
- pix_type='rectangular', apply_derotation=False)
-
- # storing the pixel mask for later use
- new_geom.mask = square_mask
-
- # create a transfer map by enumerating all pixel positions in a 2D histogram
- hex_to_rect_map = np.histogramdd([rot_y, rot_x],
- bins=(y_edges, x_edges),
- weights=np.arange(len(signal)))[0].astype(int)
- # bins that do not correspond to the original image get an entry of `-1`
- hex_to_rect_map[~square_mask] = -1
-
- if signal.ndim > 1:
- long_map = []
- for i in range(signal.shape[-1]):
- tmp_map = hex_to_rect_map + i * (np.max(hex_to_rect_map) + 1)
- tmp_map[~square_mask] = -1
- long_map.append(tmp_map)
- hex_to_rect_map = np.array(long_map)
-
- if key is not None:
- # if a key is given, store the essential objects in a buffer
- rot_buffer[key] = (geom, new_geom, hex_to_rect_map)
-
- # done `if key in rot_buffer`
-
- # create the rotated rectangular image by applying `hex_to_rect_map` to the flat,
- # extended input image
- # `input_img_ext` is the flattened input image extended by one entry that contains NaN
- # since `hex_to_rect_map` contains `-1` for "fake" pixels, it maps this extra NaN
- # value at the last array position to any bin that does not correspond to a pixel of
- # the original image
- input_img_ext = np.full(np.prod(signal.shape) + 1, np.nan)
-
- # the way the map is produced, it has the time dimension as axis=0;
- # but `signal` has it as axis=-1, so we need to roll the axes back and forth a bit.
- # if there is no time dimension, `signal` is a 1d array and `rollaxis` has no effect.
- input_img_ext[:-1] = np.rollaxis(signal, axis=-1, start=0).ravel()
-
- # now apply the transfer map
- rot_img = input_img_ext[hex_to_rect_map]
-
- # if there is a time dimension, roll the time axis back to the last position
- try:
- rot_img = np.rollaxis(rot_img, 0, 3)
- except ValueError:
- pass
-
- return new_geom, rot_img
-
-
-def convert_geometry_rect2d_back_to_hexe1d(geom, signal, key=None, add_rot=None):
- """reverts the geometry distortion performed by convert_geometry_hexe1d_to_rect_2d
- back to a hexagonal grid stored in 1D arrays
-
- Parameters
- ----------
- geom : CameraGeometry
- geometry object where pixel positions are stored in a 2D
- rectangular camera grid
- signal : ndarray
- pixel intensity stored in a 2D rectangular camera grid
- key:
- key to retrieve buffered geometry information
- add_rot:
- not used -- only here for backwards compatibility
-
- Returns
- -------
- old_geom : CameraGeometry
- the original geometry of the image
- signal : ndarray
- 1D (no timing) or 2D (with timing) array of the pmt signals
- """
-
- if key in rot_buffer:
- (old_geom, new_geom, hex_square_map) = rot_buffer[key]
- else:
- raise KeyError("key '{}' not found in the buffer".format(key)
- + " -- don't know how to undo rotation")
-
- # the output image has as many entries as there are non-negative values in the
- # transfer map (this accounts for time as well)
- unrot_img = np.zeros(np.count_nonzero(hex_square_map >= 0))
-
- # rearrange input `signal` according to the mask and map
- # (the dots in the brackets expand the mask to account for a possible time dimension)
- # `atleast_3d` ensures that there is a third axis that we can roll to the front
- # even if there is no time; if we'd use `axis=-1` instead, in cas of no time
- # dimensions, we would rotate the x and y axes, resulting in a mirrored image
- # `squeeze` reduces the added axis again in the no-time-slices cases
- unrot_img[hex_square_map[..., new_geom.mask]] = \
- np.squeeze(np.rollaxis(np.atleast_3d(signal), 2, 0))[..., new_geom.mask]
-
- # if `signal` has a third dimension, that is the time
- # and we need to roll some axes again...
- if signal.ndim == 3:
-
- # unrot_img[hex_square_map[..., new_geom.mask]] = \
- # np.rollaxis(signal, -1, 0)[..., new_geom.mask]
-
- # reshape the image so that the time is the first axis
- # and then roll the time to the back
- unrot_img = unrot_img.reshape((signal.shape[2],
- np.count_nonzero(new_geom.mask)))
- unrot_img = np.rollaxis(unrot_img, -1, 0)
- # else:
- # unrot_img[hex_square_map[new_geom.mask]] = \
- # signal[new_geom.mask]
-
-
- return old_geom, unrot_img
-
-
-convert_geometry_1d_to_2d = convert_geometry_hex1d_to_rect2d
-convert_geometry_back = convert_geometry_rect2d_back_to_hexe1d
+from .geometry_converter_astri import astri_to_2d_array, array_2d_to_astri
+from .geometry_converter_chec import chec_to_2d_array, array_2d_to_chec
diff --git a/ctapipe/image/geometry_converter_astri.py b/ctapipe/image/geometry_converter_astri.py
new file mode 100644
index 0000000000000000000000000000000000000000..bd96875ae06fa95361a13c5bdd9cea4f43e7df32
--- /dev/null
+++ b/ctapipe/image/geometry_converter_astri.py
@@ -0,0 +1,157 @@
+__author__ = "@jeremiedecock"
+
+import numpy as np
+
+
+def astri_transformation_map():
+ """Make the transformation map to turn the 1D array of an ASTRI image into a
+ rectangular 2D array. This will add new "fake" pixels that are filled with NaN value.
+ """
+
+ # By default, pixels map to the last element of input_img_ext (i.e. NaN)
+ img_map = np.full([8 * 7, 8 * 7], -1, dtype=int)
+
+ # Map values
+ img_map[0*8:1*8, 2*8:3*8] = np.arange(64).reshape([8,8])[::-1,:] + 34 * 64
+ img_map[0*8:1*8, 3*8:4*8] = np.arange(64).reshape([8,8])[::-1,:] + 35 * 64
+ img_map[0*8:1*8, 4*8:5*8] = np.arange(64).reshape([8,8])[::-1,:] + 36 * 64
+
+ img_map[1*8:2*8, 1*8:2*8] = np.arange(64).reshape([8,8])[::-1,:] + 29 * 64
+ img_map[1*8:2*8, 2*8:3*8] = np.arange(64).reshape([8,8])[::-1,:] + 30 * 64
+ img_map[1*8:2*8, 3*8:4*8] = np.arange(64).reshape([8,8])[::-1,:] + 31 * 64
+ img_map[1*8:2*8, 4*8:5*8] = np.arange(64).reshape([8,8])[::-1,:] + 32 * 64
+ img_map[1*8:2*8, 5*8:6*8] = np.arange(64).reshape([8,8])[::-1,:] + 33 * 64
+
+ img_map[2*8:3*8, 0*8:1*8] = np.arange(64).reshape([8,8])[::-1,:] + 22 * 64
+ img_map[2*8:3*8, 1*8:2*8] = np.arange(64).reshape([8,8])[::-1,:] + 23 * 64
+ img_map[2*8:3*8, 2*8:3*8] = np.arange(64).reshape([8,8])[::-1,:] + 24 * 64
+ img_map[2*8:3*8, 3*8:4*8] = np.arange(64).reshape([8,8])[::-1,:] + 25 * 64
+ img_map[2*8:3*8, 4*8:5*8] = np.arange(64).reshape([8,8])[::-1,:] + 26 * 64
+ img_map[2*8:3*8, 5*8:6*8] = np.arange(64).reshape([8,8])[::-1,:] + 27 * 64
+ img_map[2*8:3*8, 6*8:7*8] = np.arange(64).reshape([8,8])[::-1,:] + 28 * 64
+
+ img_map[3*8:4*8, 0*8:1*8] = np.arange(64).reshape([8,8])[::-1,:] + 15 * 64
+ img_map[3*8:4*8, 1*8:2*8] = np.arange(64).reshape([8,8])[::-1,:] + 16 * 64
+ img_map[3*8:4*8, 2*8:3*8] = np.arange(64).reshape([8,8])[::-1,:] + 17 * 64
+ img_map[3*8:4*8, 3*8:4*8] = np.arange(64).reshape([8,8])[::-1,:] + 18 * 64
+ img_map[3*8:4*8, 4*8:5*8] = np.arange(64).reshape([8,8])[::-1,:] + 19 * 64
+ img_map[3*8:4*8, 5*8:6*8] = np.arange(64).reshape([8,8])[::-1,:] + 20 * 64
+ img_map[3*8:4*8, 6*8:7*8] = np.arange(64).reshape([8,8])[::-1,:] + 21 * 64
+
+ img_map[4*8:5*8, 0*8:1*8] = np.arange(64).reshape([8,8])[::-1,:] + 8 * 64
+ img_map[4*8:5*8, 1*8:2*8] = np.arange(64).reshape([8,8])[::-1,:] + 9 * 64
+ img_map[4*8:5*8, 2*8:3*8] = np.arange(64).reshape([8,8])[::-1,:] + 10 * 64
+ img_map[4*8:5*8, 3*8:4*8] = np.arange(64).reshape([8,8])[::-1,:] + 11 * 64
+ img_map[4*8:5*8, 4*8:5*8] = np.arange(64).reshape([8,8])[::-1,:] + 12 * 64
+ img_map[4*8:5*8, 5*8:6*8] = np.arange(64).reshape([8,8])[::-1,:] + 13 * 64
+ img_map[4*8:5*8, 6*8:7*8] = np.arange(64).reshape([8,8])[::-1,:] + 14 * 64
+
+ img_map[5*8:6*8, 1*8:2*8] = np.arange(64).reshape([8,8])[::-1,:] + 3 * 64
+ img_map[5*8:6*8, 2*8:3*8] = np.arange(64).reshape([8,8])[::-1,:] + 4 * 64
+ img_map[5*8:6*8, 3*8:4*8] = np.arange(64).reshape([8,8])[::-1,:] + 5 * 64
+ img_map[5*8:6*8, 4*8:5*8] = np.arange(64).reshape([8,8])[::-1,:] + 6 * 64
+ img_map[5*8:6*8, 5*8:6*8] = np.arange(64).reshape([8,8])[::-1,:] + 7 * 64
+
+ img_map[6*8:7*8, 2*8:3*8] = np.arange(64).reshape([8,8])[::-1,:] + 0 * 64
+ img_map[6*8:7*8, 3*8:4*8] = np.arange(64).reshape([8,8])[::-1,:] + 1 * 64
+ img_map[6*8:7*8, 4*8:5*8] = np.arange(64).reshape([8,8])[::-1,:] + 2 * 64
+
+ return img_map
+
+
+def astri_to_2d_array_no_crop(input_img, img_map=astri_transformation_map()):
+ """Convert images comming form "ASTRI" telescopes in order to get regular 2D
+ "rectangular" images directly usable with most image processing tools.
+
+ Parameters
+ ----------
+ input_img : numpy.array
+ The image to convert
+
+ Returns
+ -------
+ A numpy.array containing the cropped image.
+ """
+
+ # Check the image
+ if len(input_img) != (37 * 64):
+ raise ValueError("The input image is not a valide ASTRI telescope image.")
+
+ # Copy the input flat ctapipe image and add one element with the NaN value in the end
+
+ input_img_ext = np.zeros(input_img.shape[0] + 1)
+ input_img_ext[:-1] = input_img[:]
+ input_img_ext[-1] = np.nan
+
+ # Make the output image
+ img_2d = input_img_ext[img_map]
+
+ return img_2d
+
+
+# for archaic consistency; there was another version of the function that would cut
+# off over-standing parts of the image
+astri_to_2d_array = astri_to_2d_array_no_crop
+
+
+def array_2d_to_astri(img_2d):
+ """Transforms the 2D ASTRI image back to 1D
+
+ Parameters
+ ----------
+ img_2d : 2d numpy array
+ the 2D ASTRI image
+
+ Returns
+ -------
+ img_1d : 1d numpy array
+ the 1D ASTRI image
+ """
+
+ img_1d = np.concatenate([
+ img_2d[6*8:7*8, 2*8:3*8][::-1,:].ravel(),
+ img_2d[6*8:7*8, 3*8:4*8][::-1,:].ravel(),
+ img_2d[6*8:7*8, 4*8:5*8][::-1,:].ravel(),
+ #
+ img_2d[5*8:6*8, 1*8:2*8][::-1,:].ravel(),
+ img_2d[5*8:6*8, 2*8:3*8][::-1,:].ravel(),
+ img_2d[5*8:6*8, 3*8:4*8][::-1,:].ravel(),
+ img_2d[5*8:6*8, 4*8:5*8][::-1,:].ravel(),
+ img_2d[5*8:6*8, 5*8:6*8][::-1,:].ravel(),
+ #
+ img_2d[4*8:5*8, 0*8:1*8][::-1,:].ravel(),
+ img_2d[4*8:5*8, 1*8:2*8][::-1,:].ravel(),
+ img_2d[4*8:5*8, 2*8:3*8][::-1,:].ravel(),
+ img_2d[4*8:5*8, 3*8:4*8][::-1,:].ravel(),
+ img_2d[4*8:5*8, 4*8:5*8][::-1,:].ravel(),
+ img_2d[4*8:5*8, 5*8:6*8][::-1,:].ravel(),
+ img_2d[4*8:5*8, 6*8:7*8][::-1,:].ravel(),
+ #
+ img_2d[3*8:4*8, 0*8:1*8][::-1,:].ravel(),
+ img_2d[3*8:4*8, 1*8:2*8][::-1,:].ravel(),
+ img_2d[3*8:4*8, 2*8:3*8][::-1,:].ravel(),
+ img_2d[3*8:4*8, 3*8:4*8][::-1,:].ravel(),
+ img_2d[3*8:4*8, 4*8:5*8][::-1,:].ravel(),
+ img_2d[3*8:4*8, 5*8:6*8][::-1,:].ravel(),
+ img_2d[3*8:4*8, 6*8:7*8][::-1,:].ravel(),
+ #
+ img_2d[2*8:3*8, 0*8:1*8][::-1,:].ravel(),
+ img_2d[2*8:3*8, 1*8:2*8][::-1,:].ravel(),
+ img_2d[2*8:3*8, 2*8:3*8][::-1,:].ravel(),
+ img_2d[2*8:3*8, 3*8:4*8][::-1,:].ravel(),
+ img_2d[2*8:3*8, 4*8:5*8][::-1,:].ravel(),
+ img_2d[2*8:3*8, 5*8:6*8][::-1,:].ravel(),
+ img_2d[2*8:3*8, 6*8:7*8][::-1,:].ravel(),
+ #
+ img_2d[1*8:2*8, 1*8:2*8][::-1,:].ravel(),
+ img_2d[1*8:2*8, 2*8:3*8][::-1,:].ravel(),
+ img_2d[1*8:2*8, 3*8:4*8][::-1,:].ravel(),
+ img_2d[1*8:2*8, 4*8:5*8][::-1,:].ravel(),
+ img_2d[1*8:2*8, 5*8:6*8][::-1,:].ravel(),
+ #
+ img_2d[0*8:1*8, 2*8:3*8][::-1,:].ravel(),
+ img_2d[0*8:1*8, 3*8:4*8][::-1,:].ravel(),
+ img_2d[0*8:1*8, 4*8:5*8][::-1,:].ravel()
+ ])
+
+ return img_1d
diff --git a/ctapipe/image/geometry_converter_chec.py b/ctapipe/image/geometry_converter_chec.py
new file mode 100644
index 0000000000000000000000000000000000000000..4219eb54650143b0cfaea3dec81af4732acb3d55
--- /dev/null
+++ b/ctapipe/image/geometry_converter_chec.py
@@ -0,0 +1,54 @@
+__author__ = "@jeremiedecock"
+
+import numpy as np
+
+
+def chec_transformation_map():
+ # By default, pixels map to the last element of input_img_ext (i.e. NaN)
+ img_map = np.full([8 * 6, 8 * 6], -1, dtype=int)
+
+ # Map values
+ img_map[:8, 8:-8] = np.arange(8 * 8 * 4).reshape([8, 8 * 4])
+ img_map[8:40, :] = np.arange(32 * 48).reshape([32, 48]) + 256
+ img_map[-8:, 8:-8] = np.arange(8 * 8 * 4).reshape([8, 8 * 4]) + 1792
+
+ return img_map
+
+
+def chec_to_2d_array(input_img, img_map=chec_transformation_map()):
+ """
+ Convert images comming form "CHEC" cameras in order to get regular 2D "rectangular"
+ images directly usable with most image processing tools.
+
+ Parameters
+ ----------
+ input_img : numpy.array
+ The image to convert
+
+ Returns
+ -------
+ A numpy.array containing the cropped image.
+ """
+
+ # Check the image
+ if len(input_img) != 2048:
+ raise ValueError("The input image is not a valide CHEC camera image.")
+
+ # Copy the input flat ctapipe image and add one element with the NaN value in the end
+
+ input_img_ext = np.zeros(input_img.shape[0] + 1)
+ input_img_ext[:-1] = input_img[:]
+ input_img_ext[-1] = np.nan
+
+ # Make the output image
+ img_2d = input_img_ext[img_map]
+
+ return img_2d
+
+
+def array_2d_to_chec(img_2d):
+
+ # Flatten image and remove NaN values
+ img_1d = img_2d[np.isfinite(img_2d)]
+
+ return img_1d
diff --git a/ctapipe/image/geometry_converter_hex.py b/ctapipe/image/geometry_converter_hex.py
new file mode 100644
index 0000000000000000000000000000000000000000..185dff70297bbd50cc913f1d237c933b2fedf936
--- /dev/null
+++ b/ctapipe/image/geometry_converter_hex.py
@@ -0,0 +1,525 @@
+__author__ = "@tino-michael"
+
+import logging
+from collections import namedtuple
+import numpy as np
+from astropy import units as u
+from numba import jit
+
+from ctapipe.instrument import CameraGeometry
+
+logger = logging.getLogger(__name__)
+
+
+__all__ = [
+ "convert_geometry_hex1d_to_rect2d",
+ "convert_geometry_rect2d_back_to_hexe1d"
+]
+
+
+def unskew_hex_pixel_grid(pix_x, pix_y, cam_angle=0 * u.deg,
+ base_angle=60 * u.deg):
+ r"""transform the pixel coordinates of a hexagonal image into an
+ orthogonal image
+
+ Parameters
+ ----------
+ pix_x, pix_y : 1D numpy arrays
+ the list of x and y coordinates of the hexagonal pixel grid
+ cam_angle : astropy.Quantity (default: 0 degrees)
+ The skewing is performed along the y-axis, therefore, one of the slanted
+ base-vectors needs to be parallel to the y-axis.
+ Some camera grids are rotated in a way that this is not the case.
+ This needs to be corrected.
+ base_angle : astropy.Quantity (default: 60 degrees)
+ the skewing angle of the hex-grid. should be 60° for regular hexagons
+
+ Returns
+ -------
+ pix_x, pix_y : 1D numpy arrays
+ the list of x and y coordinates of the slanted, orthogonal pixel grid
+
+ Notes
+ -----
+ The correction on the pixel position r can be described by a rotation R around
+ one angle and a sheer S along a certain axis:
+
+ .. math::
+ r' = S \cdot R \cdot r
+
+ .. math::
+ \begin{pmatrix}
+ x' \\
+ y'
+ \end{pmatrix}
+ =
+ \begin{pmatrix}
+ 1 & 0 \\
+ -1/\tan & 1
+ \end{pmatrix}
+ \cdot
+ \begin{pmatrix}
+ \cos & -\sin \\
+ \sin & \cos
+ \end{pmatrix}
+ \cdot
+ \begin{pmatrix}
+ x \\
+ y
+ \end{pmatrix}
+
+ .. math::
+ \begin{pmatrix}
+ x' \\
+ y'
+ \end{pmatrix}
+ =
+ \begin{pmatrix}
+ \cos & -\sin \\
+ \sin-\cos/\tan & \sin/\tan+\cos
+ \end{pmatrix}
+ \cdot
+ \begin{pmatrix}
+ x \\
+ y
+ \end{pmatrix}
+
+ """
+
+ tan_angle = np.tan(base_angle)
+
+ # If the camera-rotation angle is non-zero, create a rotation+sheering
+ # matrix for the pixel coordinates
+ if cam_angle != 0 * u.deg:
+ sin_angle = np.sin(cam_angle)
+ cos_angle = np.cos(cam_angle)
+
+ # the correction on the pixel position r can be described by a
+ # rotation R around one angle and a sheer S along a certain axis:
+ #
+ # r' = S * R * r
+ # (x') = ( 1 0) * (cos -sin) * (x) = ( cos -sin ) * (x)
+ # (y') (-1/tan 1) (sin cos) (y) (sin-cos/tan sin/tan+cos) * (y)
+ rot_mat = np.array(
+ [[cos_angle, -sin_angle],
+ [sin_angle - cos_angle / tan_angle,
+ sin_angle / tan_angle + cos_angle]])
+
+ else:
+ # if we don't rotate the camera, only perform the sheer
+ rot_mat = np.array([[1, 0], [-1 / tan_angle, 1]])
+
+ rotated = np.dot(rot_mat, [pix_x.value, pix_y.value])
+ rot_x = rotated[0] * pix_x.unit
+ rot_y = rotated[1] * pix_x.unit
+ return rot_x, rot_y
+
+
+def reskew_hex_pixel_grid(pix_x, pix_y, cam_angle=0 * u.deg,
+ base_angle=60 * u.deg):
+ r"""skews the orthogonal coordinates back to the hexagonal ones
+
+ Parameters
+ ----------
+ pix_x, pix_y : 1D numpy arrays
+ the list of x and y coordinates of the slanted, orthogonal pixel grid
+ cam_angle : astropy.Quantity (default: 0 degrees)
+ The skewing is performed along the y-axis, therefore, one of the slanted
+ base-vectors needs to be parallel to the y-axis.
+ Some camera grids are rotated in a way that this is not the case.
+ This needs to be corrected.
+ base_angle : astropy.Quantity (default: 60 degrees)
+ the skewing angle of the hex-grid. should be 60° for regular hexagons
+
+ Returns
+ -------
+ pix_x, pix_y : 1D numpy arrays
+ the list of x and y coordinates of the hexagonal pixel grid
+
+ Notes
+ -----
+ To revert the rotation, we need to find matrices S' and R' with
+ :math:`S' \cdot S = 1` and :math:`R' \cdot R = 1`,
+ so that :math:`r = R' \cdot S' \cdot S \cdot R \cdot r = R' \cdot S' \cdot r'`:
+
+ .. math::
+ \begin{pmatrix}
+ x \\
+ y
+ \end{pmatrix}
+ =
+ \begin{pmatrix}
+ \cos & \sin \\
+ -\sin & \cos
+ \end{pmatrix}
+ \cdot
+ \begin{pmatrix}
+ 1 & 0 \\
+ 1/\tan & 1
+ \end{pmatrix}
+ \cdot
+ \begin{pmatrix}
+ x' \\
+ y'
+ \end{pmatrix}
+
+ .. math::
+ \begin{pmatrix}
+ x \\
+ y
+ \end{pmatrix}
+ =
+ \begin{pmatrix}
+ \cos+\sin/\tan & \sin \\
+ \cos/\tan-\sin & \cos
+ \end{pmatrix}
+ \cdot
+ \begin{pmatrix}
+ x' \\
+ y'
+ \end{pmatrix}
+
+ """
+
+ tan_angle = np.tan(base_angle)
+
+ # If the camera-rotation angle is non-zero, create a rotation+sheering
+ # matrix for the pixel coordinates
+ if cam_angle != 0 * u.deg:
+ sin_angle = np.sin(cam_angle)
+ cos_angle = np.cos(cam_angle)
+
+ # to revert the rotation, we need to find matrices S' and R'
+ # S' * S = 1 and R' * R = 1
+ # so that
+ # r = R' * S' * S * R * r = R' * S' * r'
+ #
+ # (x) = ( cos sin) * ( 1 0) * (x') = (cos+sin/tan sin) * (x')
+ # (y) (-sin cos) (1/tan 1) (y') (cos/tan-sin cos) (y')
+
+ rot_mat = np.array(
+ [[cos_angle + sin_angle / tan_angle, sin_angle],
+ [cos_angle / tan_angle - sin_angle, cos_angle]])
+
+ else:
+ # if we don't rotate the camera, only perform the sheer
+ rot_mat = np.array([[1, 0], [1 / tan_angle, 1]])
+
+ rotated = np.dot(rot_mat, [pix_x.value, pix_y.value])
+ rot_x = rotated[0] * pix_x.unit
+ rot_y = rotated[1] * pix_x.unit
+ return rot_x, rot_y
+
+
+@jit
+def reskew_hex_pixel_from_orthogonal_edges(x_edges, y_edges, square_mask):
+ """extracts and skews the pixel coordinates from a 2D orthogonal
+ histogram (i.e. the bin-edges) and skews them into the hexagonal
+ image while selecting only the pixel that are selected by the
+ given mask
+
+ Parameters
+ ----------
+ x_edges, y_edges : 1darrays
+ the bin edges of the 2D histogram
+ square_mask : 2darray
+ mask that selects the pixels actually belonging to the camera
+
+ Returns
+ -------
+ unrot_x, unrot_y : 1darrays
+ pixel coordinated reskewed into the hexagonal camera grid
+ """
+
+ unrot_x, unrot_y = [], []
+ for i, x in enumerate((x_edges[:-1] + x_edges[1:]) / 2):
+ for j, y in enumerate((y_edges[:-1] + y_edges[1:]) / 2):
+ if square_mask[i][j]:
+ x_unrot, y_unrot = reskew_hex_pixel_grid(x, y)
+ unrot_x.append(x_unrot)
+ unrot_y.append(y_unrot)
+ return unrot_x, unrot_y
+
+
+@jit
+def get_orthogonal_grid_edges(pix_x, pix_y, scale_aspect=True):
+ """calculate the bin edges of the slanted, orthogonal pixel grid to
+ resample the pixel signals with np.histogramdd right after.
+
+ Parameters
+ ----------
+ pix_x, pix_y : 1D numpy arrays
+ the list of x and y coordinates of the slanted, orthogonal pixel grid
+ scale_aspect : boolean (default: True)
+ if True, rescales the x-coordinates to create square pixels
+ (instead of rectangular ones)
+
+ Returns
+ --------
+ x_edges, y_edges : 1D numpy arrays
+ the bin edges for the slanted, orthogonal pixel grid
+ x_scale : float
+ factor by which the x-coordinates have been scaled
+ """
+
+ # finding the size of the square patches
+ d_x = 99 * u.m
+ d_y = 99 * u.m
+ x_base = pix_x[0]
+ y_base = pix_y[0]
+ for x, y in zip(pix_x, pix_y):
+ if abs(y - y_base) < abs(x - x_base):
+ d_x = min(d_x, abs(x - x_base))
+ for x, y in zip(pix_x, pix_y):
+ if abs(y - y_base) > abs(x - x_base):
+ d_y = min(d_y, abs(y - y_base))
+
+ x_scale = 1
+ if scale_aspect:
+ x_scale = d_y / d_x
+ pix_x *= x_scale
+ d_x = d_y
+
+ # with the maximal extension of the axes and the size of the pixels,
+ # determine the number of bins in each direction
+ NBinsx = np.around(abs(max(pix_x) - min(pix_x)) / d_x) + 2
+ NBinsy = np.around(abs(max(pix_y) - min(pix_y)) / d_y) + 2
+ x_edges = np.linspace(min(pix_x).value, max(pix_x).value, NBinsx)
+ y_edges = np.linspace(min(pix_y).value, max(pix_y).value, NBinsy)
+
+ return x_edges, y_edges, x_scale
+
+
+rot_buffer = {}
+
+
+def convert_geometry_hex1d_to_rect2d(geom, signal, key=None, add_rot=0):
+ """converts the geometry object of a camera with a hexagonal grid into
+ a square grid by slanting and stretching the 1D arrays of pixel x
+ and y positions and signal intensities are converted to 2D
+ arrays. If the signal array contains a time-dimension it is
+ conserved.
+
+ Parameters
+ ----------
+ geom : CameraGeometry object
+ geometry object of hexagonal cameras
+ signal : ndarray
+ 1D (no timing) or 2D (with timing) array of the pmt signals
+ key : (default: None)
+ arbitrary key to store the transformed geometry in a buffer
+ add_rot : int/float (default: 0)
+ parameter to apply an additional rotation of `add_rot` times 60°
+
+ Returns
+ -------
+ new_geom : CameraGeometry object
+ geometry object of the slanted picture now with a rectangular
+ grid and a 2D grid for the pixel positions. contains now a 2D
+ masking array signifying which of the pixels came from the
+ original geometry and which are simply fillers from the
+ rectangular grid
+ rot_img : ndarray 2D (no timing) or 3D (with timing)
+ the rectangular signal image
+ """
+
+ if key in rot_buffer:
+
+ # if the conversion with this key was done before and stored,
+ # just read it in
+ (geom, new_geom, hex_to_rect_map) = rot_buffer[key]
+ else:
+
+ # otherwise, we have to do the conversion first now,
+ # skew all the coordinates of the original geometry
+
+ # extra_rot is the angle to get back to aligned hexagons with flat
+ # tops. Note that the pixel rotation angle brings the camera so that
+ # hexagons have a point at the top, so need to go 30deg back to
+ # make them flat
+ extra_rot = geom.pix_rotation - 30 * u.deg
+
+ # total rotation angle:
+ rot_angle = (add_rot * 60 * u.deg) - extra_rot
+
+ logger.debug("geom={}".format(geom))
+ logger.debug("rot={}, extra={}".format(rot_angle, extra_rot))
+
+ rot_x, rot_y = unskew_hex_pixel_grid(geom.pix_x, geom.pix_y,
+ cam_angle=rot_angle)
+
+ # with all the coordinate points, we can define the bin edges
+ # of a 2D histogram
+ x_edges, y_edges, x_scale = get_orthogonal_grid_edges(rot_x, rot_y)
+
+ # this histogram will introduce bins that do not correspond to
+ # any pixel from the original geometry. so we create a mask to
+ # remember the true camera pixels by simply throwing all pixel
+ # positions into numpy.histogramdd: proper pixels contain the
+ # value 1, false pixels the value 0.
+ square_mask = np.histogramdd([rot_y, rot_x],
+ bins=(y_edges, x_edges))[0].astype(bool)
+
+ # to be consistent with the pixel intensity, instead of saving
+ # only the rotated positions of the true pixels (rot_x and
+ # rot_y), create 2D arrays of all x and y positions (also the
+ # false ones).
+ grid_x, grid_y = np.meshgrid((x_edges[:-1] + x_edges[1:]) / 2.,
+ (y_edges[:-1] + y_edges[1:]) / 2.)
+
+ ids = []
+ # instead of blindly enumerating all pixels, let's instead
+ # store a list of all valid -- i.e. picked by the mask -- 2D
+ # indices
+ for i, row in enumerate(square_mask):
+ for j, val in enumerate(row):
+ if val is True:
+ ids.append((i, j))
+
+ # the area of the pixels (note that this is still a deformed
+ # image)
+ pix_area = np.ones_like(grid_x) \
+ * (x_edges[1] - x_edges[0]) * (y_edges[1] - y_edges[0])
+
+ # creating a new geometry object with the attributes we just determined
+ new_geom = CameraGeometry(
+ cam_id=geom.cam_id + "_rect",
+ pix_id=ids, # this is a list of all the valid coordinate pairs now
+ pix_x=grid_x * u.m,
+ pix_y=grid_y * u.m,
+ pix_area=pix_area * u.m ** 2,
+ neighbors=geom.neighbors,
+ pix_type='rectangular', apply_derotation=False)
+
+ # storing the pixel mask for later use
+ new_geom.mask = square_mask
+
+ # create a transfer map by enumerating all pixel positions in a 2D histogram
+ hex_to_rect_map = np.histogramdd([rot_y, rot_x],
+ bins=(y_edges, x_edges),
+ weights=np.arange(len(signal)))[0].astype(int)
+ # bins that do not correspond to the original image get an entry of `-1`
+ hex_to_rect_map[~square_mask] = -1
+
+ if signal.ndim > 1:
+ long_map = []
+ for i in range(signal.shape[-1]):
+ tmp_map = hex_to_rect_map + i * (np.max(hex_to_rect_map) + 1)
+ tmp_map[~square_mask] = -1
+ long_map.append(tmp_map)
+ hex_to_rect_map = np.array(long_map)
+
+ if key is not None:
+ # if a key is given, store the essential objects in a buffer
+ rot_buffer[key] = (geom, new_geom, hex_to_rect_map)
+
+ # done `if key in rot_buffer`
+
+ # create the rotated rectangular image by applying `hex_to_rect_map` to the flat,
+ # extended input image
+ # `input_img_ext` is the flattened input image extended by one entry that contains NaN
+ # since `hex_to_rect_map` contains `-1` for "fake" pixels, it maps this extra NaN
+ # value at the last array position to any bin that does not correspond to a pixel of
+ # the original image
+ input_img_ext = np.full(np.prod(signal.shape) + 1, np.nan)
+
+ # the way the map is produced, it has the time dimension as axis=0;
+ # but `signal` has it as axis=-1, so we need to roll the axes back and forth a bit.
+ # if there is no time dimension, `signal` is a 1d array and `rollaxis` has no effect.
+ input_img_ext[:-1] = np.rollaxis(signal, axis=-1, start=0).ravel()
+
+ # now apply the transfer map
+ rot_img = input_img_ext[hex_to_rect_map]
+
+ # if there is a time dimension, roll the time axis back to the last position
+ try:
+ rot_img = np.rollaxis(rot_img, 0, 3)
+ except ValueError:
+ pass
+
+ return new_geom, rot_img
+
+
+def convert_geometry_rect2d_back_to_hexe1d(geom, signal, key=None, add_rot=None):
+ """reverts the geometry distortion performed by convert_geometry_hexe1d_to_rect_2d
+ back to a hexagonal grid stored in 1D arrays
+
+ Parameters
+ ----------
+ geom : CameraGeometry
+ geometry object where pixel positions are stored in a 2D
+ rectangular camera grid
+ signal : ndarray
+ pixel intensity stored in a 2D rectangular camera grid
+ key:
+ key to retrieve buffered geometry information
+ add_rot:
+ not used -- only here for backwards compatibility
+
+ Returns
+ -------
+ old_geom : CameraGeometry
+ the original geometry of the image
+ signal : ndarray
+ 1D (no timing) or 2D (with timing) array of the pmt signals
+
+ Note
+ ----
+ The back-conversion works with an internal buffer to store the transfer map (which
+ was produced in the first conversion). If `key` is not found in said buffer, this
+ function tries to perform a mock conversion. For this, it needs a `CameraGeometry`
+ instance of the original camera layout, which it tries to load by name (i.e.
+ the `cam_id`). The function assumes the original `cam_id` can be inferred from the
+ given, modified one by: `geom.cam_id.split('_')[0]`.
+ """
+
+
+ if key not in rot_buffer:
+ # if the key is not in the buffer from the initial conversion (maybe because you
+ # did it in another process?), perform a mock conversion here
+ # ATTENTION assumes the original cam_id can be inferred from the given, modified
+ # one by by `geom.cam_id.split('_')[0]`
+ try:
+ orig_geom = CameraGeometry.from_name(geom.cam_id.split('_')[0])
+ except:
+ raise ValueError("could not deduce `CameraGeometry` from given `geom`...\n"
+ "please provide a `geom`, so that "
+ "`geom.cam_id.split('_')[0]` is a known `cam_id`")
+
+ orig_signal = np.zeros(len(orig_geom.pix_x))
+ convert_geometry_hex1d_to_rect2d(geom=orig_geom, signal=orig_signal,
+ key=key, add_rot=add_rot)
+
+ (old_geom, new_geom, hex_square_map) = rot_buffer[key]
+
+ # the output image has as many entries as there are non-negative values in the
+ # transfer map (this accounts for time as well)
+ unrot_img = np.zeros(np.count_nonzero(hex_square_map >= 0))
+
+ # rearrange input `signal` according to the mask and map
+ # (the dots in the brackets expand the mask to account for a possible time dimension)
+ # `atleast_3d` ensures that there is a third axis that we can roll to the front
+ # even if there is no time; if we'd use `axis=-1` instead, in cas of no time
+ # dimensions, we would rotate the x and y axes, resulting in a mirrored image
+ # `squeeze` reduces the added axis again in the no-time-slices cases
+ unrot_img[hex_square_map[..., new_geom.mask]] = \
+ np.squeeze(np.rollaxis(np.atleast_3d(signal), 2, 0))[..., new_geom.mask]
+
+ # if `signal` has a third dimension, that is the time
+ # and we need to roll some axes again...
+ if signal.ndim == 3:
+
+ # unrot_img[hex_square_map[..., new_geom.mask]] = \
+ # np.rollaxis(signal, -1, 0)[..., new_geom.mask]
+
+ # reshape the image so that the time is the first axis
+ # and then roll the time to the back
+ unrot_img = unrot_img.reshape((signal.shape[2],
+ np.count_nonzero(new_geom.mask)))
+ unrot_img = np.rollaxis(unrot_img, -1, 0)
+ # else:
+ # unrot_img[hex_square_map[new_geom.mask]] = \
+ # signal[new_geom.mask]
+
+
+ return old_geom, unrot_img
diff --git a/ctapipe/image/tests/test_geometry_converter.py b/ctapipe/image/tests/test_geometry_converter.py
index 70c540d36cc85cce2acf1ab4b57a4c72f64dbe9b..87b1ba643b0d7d2611940532bdeb7c38260db7f3 100644
--- a/ctapipe/image/tests/test_geometry_converter.py
+++ b/ctapipe/image/tests/test_geometry_converter.py
@@ -4,8 +4,10 @@ import numpy as np
from matplotlib import pyplot as plt
from ctapipe.image.cleaning import tailcuts_clean
-from ctapipe.image.geometry_converter import convert_geometry_1d_to_2d, \
- convert_geometry_back
+from ctapipe.image.geometry_converter import (convert_geometry_hex1d_to_rect2d,
+ convert_geometry_rect2d_back_to_hexe1d,
+ astri_to_2d_array, array_2d_to_astri,
+ chec_to_2d_array, array_2d_to_chec)
from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError
from ctapipe.instrument import CameraGeometry
from ctapipe.io.hessio import hessio_event_source
@@ -24,9 +26,6 @@ def test_convert_geometry(cam_id, rot):
geom = CameraGeometry.from_name(cam_id)
- if geom.pix_type == 'rectangular':
- return # skip non-hexagonal cameras, since they don't need conversion
-
model = generate_2d_shower_model(centroid=(0.4, 0), width=0.01, length=0.03,
psi="25d")
@@ -36,12 +35,30 @@ def test_convert_geometry(cam_id, rot):
hillas_0 = hillas_parameters(geom, image)
- geom2d, image2d = convert_geometry_1d_to_2d(geom, image,
+ if geom.pix_type == 'hexagonal':
+ convert_geometry_1d_to_2d = convert_geometry_hex1d_to_rect2d
+ convert_geometry_back = convert_geometry_rect2d_back_to_hexe1d
+
+ geom2d, image2d = convert_geometry_1d_to_2d(geom, image,
+ geom.cam_id + str(rot),
+ add_rot=rot)
+ geom1d, image1d = convert_geometry_back(geom2d, image2d,
geom.cam_id + str(rot),
add_rot=rot)
- geom1d, image1d = convert_geometry_back(geom2d, image2d,
- geom.cam_id + str(rot),
- add_rot=rot)
+
+ else:
+ if geom.cam_id == "ASTRICam":
+ convert_geometry_1d_to_2d = astri_to_2d_array
+ convert_geometry_back = array_2d_to_astri
+ elif geom.cam_id == "CHEC":
+ convert_geometry_1d_to_2d = chec_to_2d_array
+ convert_geometry_back = array_2d_to_chec
+ else:
+ print("camera {geom.cam_id} not implemented")
+ return
+
+ image2d = convert_geometry_1d_to_2d(image)
+ image1d = convert_geometry_back(image2d)
hillas_1 = hillas_parameters(geom, image1d)
@@ -54,6 +71,42 @@ def test_convert_geometry(cam_id, rot):
# TODO: test other parameters
+@pytest.mark.parametrize("rot", [3, ])
+@pytest.mark.parametrize("cam_id", cam_ids)
+def test_convert_geometry_mock(cam_id, rot):
+ """here we use a different key for the back conversion to trigger the mock conversion
+ """
+
+ geom = CameraGeometry.from_name(cam_id)
+
+ model = generate_2d_shower_model(centroid=(0.4, 0), width=0.01, length=0.03,
+ psi="25d")
+
+ _, image, _ = make_toymodel_shower_image(geom, model.pdf,
+ intensity=50,
+ nsb_level_pe=100)
+
+ hillas_0 = hillas_parameters(geom, image)
+
+ if geom.pix_type == 'hexagonal':
+ convert_geometry_1d_to_2d = convert_geometry_hex1d_to_rect2d
+ convert_geometry_back = convert_geometry_rect2d_back_to_hexe1d
+
+ geom2d, image2d = convert_geometry_1d_to_2d(geom, image, key=None,
+ add_rot=rot)
+ geom1d, image1d = convert_geometry_back(geom2d, image2d,
+ "_".join([geom.cam_id,
+ str(rot), "mock"]),
+ add_rot=rot)
+ else:
+ # originally rectangular geometries don't need a buffer and therefore no mock
+ # conversion
+ return
+
+ hillas_1 = hillas_parameters(geom, image1d)
+ assert np.abs(hillas_1.phi - hillas_0.phi).deg < 1.0
+
+
# def plot_cam(geom, geom2d, geom1d, image, image2d, image1d):
# # plt.viridis()
# plt.figure(figsize=(12, 4))
diff --git a/ctapipe/instrument/tests/test_camera.py b/ctapipe/instrument/tests/test_camera.py
index faae48c10d1c03c000e4959fe20c9eda98a61502..0f32d88038accf6f3a622cf9a9cea4e18c765783 100644
--- a/ctapipe/instrument/tests/test_camera.py
+++ b/ctapipe/instrument/tests/test_camera.py
@@ -3,25 +3,22 @@ from astropy import units as u
from ctapipe.instrument import CameraGeometry
from ctapipe.instrument.camera import _find_neighbor_pixels, \
_get_min_pixel_seperation
-from numpy import median
import pytest
cam_ids = CameraGeometry.get_known_camera_names()
-def test_load_by_name():
+def test_known_camera_names():
cams = CameraGeometry.get_known_camera_names()
assert len(cams) > 4
assert 'FlashCam' in cams
assert 'NectarCam' in cams
-
for cam in cams:
geom = CameraGeometry.from_name(cam)
geom.info()
-
def test_make_rectangular_camera_geometry():
geom = CameraGeometry.make_rectangular()
assert(geom.pix_x.shape == geom.pix_y.shape)
@@ -35,7 +32,7 @@ def test_load_hess_camera():
def test_guess_camera():
px = np.linspace(-10, 10, 11328) * u.m
py = np.linspace(-10, 10, 11328) * u.m
- geom = CameraGeometry.guess(px, py,0 * u.m)
+ geom = CameraGeometry.guess(px, py, 0 * u.m)
assert geom.pix_type.startswith('rect')
@@ -75,7 +72,6 @@ def test_neighbor_pixels(cam_id):
assert n_neighbors.count(1) == 0 # no pixel should have a single neighbor
-
def test_to_and_from_table():
geom = CameraGeometry.from_name("LSTCam")
tab = geom.to_table()
@@ -102,11 +98,6 @@ def test_write_read(tmpdir):
assert (geom.pix_area == geom2.pix_area).all()
assert geom.pix_type == geom2.pix_type
-def test_known_cameras():
- cams = CameraGeometry.get_known_camera_names()
- assert 'FlashCam' in cams
- assert len(cams) > 3
-
def test_precal_neighbors():
geom = CameraGeometry(cam_id="TestCam",
@@ -114,10 +105,12 @@ def test_precal_neighbors():
pix_x=np.arange(3)*u.deg,
pix_y=np.arange(3)*u.deg,
pix_area=np.ones(3)*u.deg**2,
- neighbors=[[1,],[0,2],[1,]],
+ neighbors=[
+ [1, ], [0, 2], [1, ]
+ ],
pix_type='rectangular',
pix_rotation="0deg",
- cam_rotation="0deg" )
+ cam_rotation="0deg")
neigh = geom.neighbors
assert len(neigh) == len(geom.pix_x)
@@ -126,7 +119,6 @@ def test_precal_neighbors():
assert nmat.shape == (len(geom.pix_x), len(geom.pix_x))
-
def test_slicing():
geom = CameraGeometry.from_name("NectarCam")
sliced1 = geom[100:200]
@@ -136,7 +128,7 @@ def test_slicing():
assert len(sliced1.pix_area) == 100
assert len(sliced1.pix_id) == 100
- sliced2 = geom[[5,7,8,9,10]]
+ sliced2 = geom[[5, 7, 8, 9, 10]]
assert sliced2.pix_id[0] == 5
assert sliced2.pix_id[1] == 7
assert len(sliced2.pix_x) == 5
diff --git a/ctapipe/instrument/tests/test_optics.py b/ctapipe/instrument/tests/test_optics.py
index 377997ffe8add1c810ce75101ab4bca192ca344b..8e090a0639694c13e5660c13463861e811337eac 100644
--- a/ctapipe/instrument/tests/test_optics.py
+++ b/ctapipe/instrument/tests/test_optics.py
@@ -1,7 +1,8 @@
-from ctapipe.instrument.optics import OpticsDescription
+from ctapipe.instrument.optics import OpticsDescription
from astropy import units as u
import pytest
+
def test_guess_optics():
od = OpticsDescription.guess(28.0*u.m)
@@ -12,6 +13,4 @@ def test_guess_optics():
assert od.mirror_type == 'DC'
with pytest.raises(KeyError):
- od2 = OpticsDescription.guess(0*u.m)
-
-
+ OpticsDescription.guess(0*u.m)
diff --git a/ctapipe/instrument/tests/test_subarray.py b/ctapipe/instrument/tests/test_subarray.py
index 774510f676fd59285e9ae78eaa5b843ad4694a1a..f8eb5cc56949efabc861ad04d3bcc44eef7c5d54 100644
--- a/ctapipe/instrument/tests/test_subarray.py
+++ b/ctapipe/instrument/tests/test_subarray.py
@@ -1,7 +1,11 @@
-from ctapipe.instrument import SubarrayDescription, TelescopeDescription, CameraGeometry
+from ctapipe.instrument import (
+ SubarrayDescription,
+ TelescopeDescription,
+)
import numpy as np
from astropy import units as u
+
def test_subarray_description():
pos = {}
@@ -11,11 +15,8 @@ def test_subarray_description():
pix_y = np.arange(1764, dtype=np.float) * u.m
for ii in range(10):
-
tel[ii] = TelescopeDescription.guess(pix_x, pix_y, foclen)
- pos[ii] = (np.random.uniform(200, size=2)-100) * u.m
-
-
+ pos[ii] = np.random.uniform(-100, 100, size=2) * u.m
sub = SubarrayDescription("test array",
tel_positions=pos,
@@ -27,6 +28,6 @@ def test_subarray_description():
assert sub.tel[0].camera is not None
assert len(sub.to_table()) == 10
- subsub = sub.select_subarray("newsub", [1,2,3,4])
+ subsub = sub.select_subarray("newsub", [1, 2, 3, 4])
assert subsub.num_tels == 4
- assert set(subsub.tels.keys()) == {1,2,3,4}
\ No newline at end of file
+ assert set(subsub.tels.keys()) == {1, 2, 3, 4}
diff --git a/ctapipe/instrument/tests/test_telescope.py b/ctapipe/instrument/tests/test_telescope.py
index d6bdbc934e5926124626fc6575520ab252f5802e..ecb003fcfa754aa0c441f17a7024aaadebf56c52 100644
--- a/ctapipe/instrument/tests/test_telescope.py
+++ b/ctapipe/instrument/tests/test_telescope.py
@@ -2,6 +2,7 @@ from ctapipe.instrument import TelescopeDescription
from astropy import units as u
import numpy as np
+
def test_telescope_description():
# setup a dummy telescope that look like an MST with FlashCam
@@ -13,4 +14,4 @@ def test_telescope_description():
assert tel.camera.cam_id == 'FlashCam'
assert tel.optics.tel_type == 'MST'
- assert str(tel) == 'MST:FlashCam'
\ No newline at end of file
+ assert str(tel) == 'MST:FlashCam'
diff --git a/ctapipe/reco/HillasReconstructor.py b/ctapipe/reco/HillasReconstructor.py
index ab6272cf0d2fcc7fd4d4df04e2141d67269a147d..d8623355892b4b362f6048f02adceb1473c6eea3 100644
--- a/ctapipe/reco/HillasReconstructor.py
+++ b/ctapipe/reco/HillasReconstructor.py
@@ -1,15 +1,15 @@
-""""
+"""
Line-intersection-based fitting.
Contact: Tino Michael <Tino.Michael@cea.fr>
"""
+
from ctapipe.utils import linalg
+from ctapipe.coordinates import coordinate_transformations as trafo
from ctapipe.reco.reco_algorithms import Reconstructor
from ctapipe.io.containers import ReconstructedShowerContainer
-from astropy.utils.decorators import deprecated
-
from itertools import combinations
import numpy as np
@@ -21,66 +21,14 @@ u.dimless = u.dimensionless_unscaled
__all__ = ['HillasReconstructor',
- 'TooFewTelescopesException',
+ 'TooFewTelescopes',
'dist_to_traces', 'MEst', 'GreatCircle']
-class TooFewTelescopesException(Exception):
+class TooFewTelescopes(Exception):
pass
-@deprecated(0.1, "will be replaced with real coord transform")
-def guess_pix_direction(pix_x, pix_y, tel_phi, tel_theta, tel_foclen):
- """
- TODO replace with proper implementation
- calculates the direction vector of corresponding to a
- (x,y) position on the camera
-
- beta is the pixel's angular distance to the centre
- according to beta / tel_view = r / maxR
- alpha is the polar angle between the y-axis and the pixel
- to find the direction the pixel is looking at:
-
- - the pixel direction is set to the telescope direction
- - offset by beta towards up
- - rotated around the telescope direction by the angle alpha
-
-
- Parameters
- -----------
- pix_x, pix_y : ndarray
- lists of x and y positions on the camera
- tel_phi, tel_theta: astropy quantities
- two angles that describe the orientation of the telescope
- tel_foclen : astropy quantity
- focal length of the telescope
-
- Returns
- -------
- pix_dirs : ndarray
- shape (n,3) list of "direction vectors"
- corresponding to a position on the camera
-
- """
-
- pix_alpha = np.arctan2(pix_y, pix_x)
-
- pix_rho = (pix_x ** 2 + pix_y ** 2) ** .5
-
- pix_beta = pix_rho / tel_foclen * u.rad
-
- tel_dir = linalg.set_phi_theta(tel_phi, tel_theta)
-
- pix_dirs = []
- for a, b in zip(pix_alpha, pix_beta):
- pix_dir = linalg.set_phi_theta(tel_phi, tel_theta - b)
-
- pix_dir = linalg.rotate_around_axis(pix_dir, tel_dir, 90 * u.deg - a)
- pix_dirs.append(pix_dir * u.dimless)
-
- return pix_dirs
-
-
def dist_to_traces(core, circles):
"""This function calculates the M-Estimator from the distances of the
suggested core position to all traces of the given GreatCircles.
@@ -244,8 +192,8 @@ class HillasReconstructor(Reconstructor):
result = ReconstructedShowerContainer()
phi, theta = linalg.get_phi_theta(dir).to(u.deg)
- # TODO make sure az and phi turn in same direction...
- result.alt, result.az = 90 * u.deg - theta, 90 * u.deg - phi
+ # TODO fix coordinates!
+ result.alt, result.az = 90 * u.deg - theta, -phi
result.core_x = pos[0]
result.core_y = pos[1]
result.core_uncert = err_est_pos
@@ -288,7 +236,7 @@ class HillasReconstructor(Reconstructor):
foclen = subarray.tel[tel_id].optics.equivalent_focal_length
circle = GreatCircle(
- guess_pix_direction(
+ trafo.pixel_position_to_direction(
np.array([moments.cen_x / u.m, p2_x / u.m]) * u.m,
np.array([moments.cen_y / u.m, p2_y / u.m]) * u.m,
tel_phi[tel_id], tel_theta[tel_id], foclen),
@@ -328,6 +276,8 @@ class HillasReconstructor(Reconstructor):
result = linalg.normalise(np.sum(crossings, axis=0)) * u.dimless
off_angles = [linalg.angle(result, cross) / u.rad for cross in crossings]
err_est_dir = np.mean(off_angles) * u.rad
+ err_est_dir = np.average(
+ off_angles, weights=[len(cross) for cross in crossings]) * u.rad
# averaging over the solutions of all permutations
return result, err_est_dir
@@ -531,7 +481,7 @@ class HillasReconstructor(Reconstructor):
dirs = []
for tel_id, hillas in hillas_dict.items():
foclen = subarray.tel[tel_id].optics.equivalent_focal_length
- max_dir, = guess_pix_direction(
+ max_dir, = trafo.pixel_position_to_direction(
np.array([hillas.cen_x / u.m]) * u.m,
np.array([hillas.cen_y / u.m]) * u.m,
tel_phi[tel_id], tel_theta[tel_id], foclen)
diff --git a/ctapipe/utils/__init__.py b/ctapipe/utils/__init__.py
index 2dc752c08665eda4153367f4a6bfd43fe8d7d72c..f5daf86c20a752a5979f260a0377ed231206b4b6 100644
--- a/ctapipe/utils/__init__.py
+++ b/ctapipe/utils/__init__.py
@@ -3,4 +3,5 @@ from .fitshistogram import Histogram
from .json2fits import json_to_fits
from .dynamic_class import dynamic_class_from_module
from .table_interpolator import TableInterpolator
-from .datasets import find_all_matching_datasets, get_table_dataset, get_dataset, find_in_path
+from .datasets import (find_all_matching_datasets, get_table_dataset, get_dataset,
+ find_in_path)
diff --git a/ctapipe/utils/linalg.py b/ctapipe/utils/linalg.py
index 9232b430786d0cf5773e44c11ad6c9fe0c96f843..1afb1ecc8490a37447a217be2c01a198a4f4a750 100644
--- a/ctapipe/utils/linalg.py
+++ b/ctapipe/utils/linalg.py
@@ -20,7 +20,7 @@ def rotation_matrix_2d(angle):
def rotate_around_axis(vec, axis, angle):
""" rotates a vector around an axis by an angle
- creates a rotation matrix and multiplies
+ creates a rotation matrix and multiplies
the initial vector with it
Parameters
@@ -107,8 +107,10 @@ def set_phi_theta_r(phi, theta, r=1):
return np.array([sin(theta) * cos(phi),
sin(theta) * sin(phi),
cos(theta)]) * r
-""" simple alias for set_phi_theta_r with default (unitless) r argument """
-set_phi_theta = lambda x, y: set_phi_theta_r(x, y)
+
+
+""" simple alias for set_phi_theta_r """
+set_phi_theta = set_phi_theta_r
def get_phi_theta(vec):
diff --git a/docs/calib/index_camera.rst b/docs/calib/index_camera.rst
index c4bb07f132a292f2171acdfb7995d716c2f35387..a3ace2be651352f343ead4e1d40992432747c924 100644
--- a/docs/calib/index_camera.rst
+++ b/docs/calib/index_camera.rst
@@ -24,7 +24,8 @@ skipped and the target data level container is unchanged (for example
performing the calibration in r1.py on data read into the R1 container will
leave the container unchanged, as there is no data in the R0 container)
-See the CTA Data Level Models definition document for information about the
+See the `CTA High-Level Data Model Definitions SYS-QA/160517
+<https://jama.cta-observatory.org/perspective.req?projectId=6&docId=26528>`_ document (CTA internal) for information about the
different data levels.
Reference/API
diff --git a/environment.yml b/environment.yml
index c0f5570284998f318ef95fcdc9de64ad41e80b6f..97085b7b2b93c12630299ac60034482d4b5de613 100644
--- a/environment.yml
+++ b/environment.yml
@@ -1,4 +1,4 @@
-name: ctadev
+name: cta-dev
channels:
- cta-observatory
- conda-forge