Skip to content

Commit

Permalink
Fixes for astropy 4.x (#1212)
Browse files Browse the repository at this point in the history 
* fix parts of geometry_converter that relied on astropy units

in astropy 4, and the latest NumPy, numpy functions *retain* units correctly
whereas before they were stripped off. Functions that relied on them disappearing
then break.

Also first fixed, and then, removed the @jit (originally it did not
ever compile in no_python mode), however now it seems it introduces
some rounding error that breaks the hex converter (probably means the
algorithm is unstable, but that's for another day!)

* a few unit fixes for ArrayDisplay

* fixed some more unit problems in ArrayDisplay

and reformatted

* reformatting
  • Loading branch information
@kosack @maxnoe
kosack authored and maxnoe committed Jan 17, 2020
1 parent 1e62e68 commit 01c715d
Show file tree
Hide file tree
Showing 5 changed files with 157 additions and 127 deletions.
110 changes: 60 additions & 50 deletions ctapipe/image/geometry_converter_hex.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,20 +4,15 @@

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"
]
__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):
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
Expand Down Expand Up @@ -100,9 +95,11 @@ def unskew_hex_pixel_grid(pix_x, pix_y, cam_angle=0 * u.deg,
# (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]])
[
[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
Expand All @@ -114,8 +111,7 @@ def unskew_hex_pixel_grid(pix_x, pix_y, cam_angle=0 * u.deg,
return rot_x, rot_y


def reskew_hex_pixel_grid(pix_x, pix_y, cam_angle=0 * u.deg,
base_angle=60 * u.deg):
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
Expand Down Expand Up @@ -197,8 +193,11 @@ def reskew_hex_pixel_grid(pix_x, pix_y, cam_angle=0 * u.deg,
# (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]])
[
[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
Expand All @@ -210,7 +209,6 @@ def reskew_hex_pixel_grid(pix_x, pix_y, cam_angle=0 * u.deg,
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
Expand Down Expand Up @@ -240,7 +238,6 @@ def reskew_hex_pixel_from_orthogonal_edges(x_edges, y_edges, square_mask):
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.
Expand All @@ -249,6 +246,7 @@ def get_orthogonal_grid_edges(pix_x, pix_y, scale_aspect=True):
----------
pix_x, pix_y : 1D numpy arrays
the list of x and y coordinates of the slanted, orthogonal pixel grid
units should be in meters, and stripped off
scale_aspect : boolean (default: True)
if True, rescales the x-coordinates to create square pixels
(instead of rectangular ones)
Expand All @@ -262,8 +260,9 @@ def get_orthogonal_grid_edges(pix_x, pix_y, scale_aspect=True):
"""

# finding the size of the square patches
d_x = 99 * u.meter # TODO: @jit may have troubles interpreting astropy.Quantities
d_y = 99 * u.meter

d_x = 99
d_y = 99
x_base = pix_x[0]
y_base = pix_y[0]
for x, y in zip(pix_x, pix_y):
Expand All @@ -284,12 +283,12 @@ def get_orthogonal_grid_edges(pix_x, pix_y, scale_aspect=True):

# with the maximal extension of the axes and the size of the pixels,
# determine the number of bins in each direction
n_bins_x = (np.around(abs(max(pix_x) - min(pix_x)) / d_x) + 2).astype(int)
n_bins_y = (np.around(abs(max(pix_y) - min(pix_y)) / d_y) + 2).astype(int)
x_edges = np.linspace(min(pix_x).value, max(pix_x).value, n_bins_x)
y_edges = np.linspace(min(pix_y).value, max(pix_y).value, n_bins_y)
n_bins_x = int(np.round_(np.abs(np.max(pix_x) - np.min(pix_x)) / d_x) + 2)
n_bins_y = int(np.round_(np.abs(np.max(pix_y) - np.min(pix_y)) / d_y) + 2)
x_edges = np.linspace(pix_x.min(), pix_x.max(), n_bins_x)
y_edges = np.linspace(pix_y.min(), pix_y.max(), n_bins_y)

return x_edges, y_edges, x_scale
return (x_edges, y_edges, x_scale)


rot_buffer = {}
Expand Down Expand Up @@ -356,27 +355,32 @@ def convert_geometry_hex1d_to_rect2d(geom, signal, key=None, add_rot=0):
logger.debug(f"geom={geom}")
logger.debug(f"rot={rot_angle}, extra={extra_rot}")

rot_x, rot_y = unskew_hex_pixel_grid(geom.pix_x, geom.pix_y,
cam_angle=rot_angle)
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)
x_edges, y_edges, x_scale = get_orthogonal_grid_edges(
rot_x.to_value(u.m), rot_y.to_value(u.m)
)

# 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)
square_mask = np.histogramdd(
[rot_y.to_value(u.m), rot_x.to_value(u.m)], 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.)
grid_x, grid_y = np.meshgrid(
(x_edges[:-1] + x_edges[1:]) / 2.0, (y_edges[:-1] + y_edges[1:]) / 2.0
)

ids = []
# instead of blindly enumerating all pixels, let's instead
Expand All @@ -389,29 +393,32 @@ def convert_geometry_hex1d_to_rect2d(geom, signal, key=None, add_rot=0):

# 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]))
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=u.Quantity(grid_x.ravel(), u.meter),
pix_y=u.Quantity(grid_y.ravel(), u.meter),
pix_x=u.Quantity(grid_x.ravel(), u.meter),
pix_y=u.Quantity(grid_y.ravel(), u.meter),
pix_area=pix_area * u.meter ** 2,
neighbors=geom.neighbors,
sampling_rate=geom.sampling_rate,
pix_type='rectangular', apply_derotation=False)
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)
hex_to_rect_map = np.histogramdd(
[rot_y.to_value(u.m), rot_x.to_value(u.m)],
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

Expand Down Expand Up @@ -454,8 +461,7 @@ def convert_geometry_hex1d_to_rect2d(geom, signal, key=None, add_rot=0):
return new_geom, rot_img


def convert_geometry_rect2d_back_to_hexe1d(geom, signal, key=None,
add_rot=None):
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
Expand Down Expand Up @@ -504,16 +510,18 @@ def convert_geometry_rect2d_back_to_hexe1d(geom, signal, key=None,
# 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])
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`")
"`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)
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]

Expand All @@ -527,8 +535,9 @@ def convert_geometry_rect2d_back_to_hexe1d(geom, signal, key=None,
# 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]
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...
Expand All @@ -538,8 +547,9 @@ def convert_geometry_rect2d_back_to_hexe1d(geom, signal, key=None,

# 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 = 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]] = \
Expand Down
54 changes: 27 additions & 27 deletions ctapipe/image/tests/test_geometry_converter.py
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,13 @@
import pytest
import numpy as np
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.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
from ctapipe.instrument import CameraGeometry
from ctapipe.image.toymodel import Gaussian
Expand All @@ -14,45 +18,43 @@


def create_mock_image(geom):
'''
"""
creates a mock image, which parameters are adapted to the camera size
'''
"""

camera_r = np.max(np.sqrt(geom.pix_x**2 + geom.pix_y**2))
camera_r = np.max(np.sqrt(geom.pix_x ** 2 + geom.pix_y ** 2))
model = Gaussian(
x=0.3 * camera_r,
y=0 * u.m,
width=0.03 * camera_r,
length=0.10 * camera_r,
psi="25d"
psi="25d",
)

_, image, _ = model.generate_image(
geom,
intensity=0.5 * geom.n_pixels,
nsb_level_pe=3,
geom, intensity=0.5 * geom.n_pixels, nsb_level_pe=3,
)
return image


@pytest.mark.parametrize("rot", [3, ])
@pytest.mark.parametrize("rot", [3,])
@pytest.mark.parametrize("cam_id", cam_ids)
def test_convert_geometry(cam_id, rot):

geom = CameraGeometry.from_name(cam_id)
image = create_mock_image(geom)
hillas_0 = hillas_parameters(geom, image)

if geom.pix_type == 'hexagonal':
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)
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
)

else:
if geom.cam_id == "ASTRICam":
Expand All @@ -79,7 +81,7 @@ def test_convert_geometry(cam_id, rot):
# TODO: test other parameters


@pytest.mark.parametrize("rot", [3, ])
@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
Expand All @@ -89,16 +91,14 @@ def test_convert_geometry_mock(cam_id, rot):
image = create_mock_image(geom)
hillas_0 = hillas_parameters(geom, image)

if geom.pix_type == 'hexagonal':
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)
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
Expand Down
Loading

0 comments on commit 01c715d

Please sign in to comment.