Skip to content

Commit

Permalink
fix extractor & update tests
Browse files Browse the repository at this point in the history
  • Loading branch information
StFroese committed Sep 15, 2023
1 parent 214213f commit 87073c6
Show file tree
Hide file tree
Showing 5 changed files with 61 additions and 187 deletions.
6 changes: 5 additions & 1 deletion ctapipe/image/extractor.py
Original file line number Diff line number Diff line change
Expand Up @@ -1126,7 +1126,11 @@ def _apply_second_pass(

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

# get the predicted times as a linear relation
Expand Down
15 changes: 8 additions & 7 deletions ctapipe/image/tests/test_concentration.py
Original file line number Diff line number Diff line change
@@ -1,14 +1,15 @@
import astropy.units as u
import pytest

from ctapipe.coordinates import TelescopeFrame
from ctapipe.image.concentration import concentration_parameters
from ctapipe.image.hillas import hillas_parameters
from ctapipe.image.tests.test_hillas import create_sample_image


def test_concentration(prod5_lst):

geom = prod5_lst.camera.geometry
geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame())
image, clean_mask = create_sample_image(psi="30d", geometry=geom)

hillas = hillas_parameters(geom[clean_mask], image[clean_mask])
Expand All @@ -22,26 +23,26 @@ def test_concentration(prod5_lst):

@pytest.mark.filterwarnings("error")
def test_width_0(prod5_lst):
geom = prod5_lst.camera.geometry
geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame())
image, clean_mask = create_sample_image(psi="30d", geometry=geom)

hillas = hillas_parameters(geom[clean_mask], image[clean_mask])
hillas.width = 0 * u.m
hillas.width = 0 * u.deg

conc = concentration_parameters(geom, image, hillas)
assert conc.core == 0


def test_no_pixels_near_cog(prod5_lst):
geom = prod5_lst.camera.geometry
geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame())
image, clean_mask = create_sample_image(psi="30d", geometry=geom)

hillas = hillas_parameters(geom[clean_mask], image[clean_mask])

# remove pixels close to cog from the cleaning pixels
clean_mask &= ((geom.pix_x - hillas.x) ** 2 + (geom.pix_y - hillas.y) ** 2) >= (
2 * geom.pixel_width**2
)
clean_mask &= (
(geom.pix_x - hillas.fov_lon) ** 2 + (geom.pix_y - hillas.fov_lat) ** 2
) >= (2 * geom.pixel_width**2)

conc = concentration_parameters(geom[clean_mask], image[clean_mask], hillas)
assert conc.cog == 0
156 changes: 36 additions & 120 deletions ctapipe/image/tests/test_hillas.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,11 @@
import numpy as np
import pytest
from astropy import units as u
from astropy.coordinates import Angle, SkyCoord
from astropy.coordinates import Angle
from numpy import isclose
from pytest import approx

from ctapipe.containers import (
CameraHillasParametersContainer,
HillasParametersContainer,
)
from ctapipe.containers import HillasParametersContainer
from ctapipe.coordinates import TelescopeFrame
from ctapipe.image import tailcuts_clean, toymodel
from ctapipe.image.hillas import HillasParameterizationError, hillas_parameters
Expand All @@ -19,17 +16,17 @@

def create_sample_image(
psi="-30d",
x=0.2 * u.m,
y=0.3 * u.m,
width=0.05 * u.m,
length=0.15 * u.m,
x=1 * u.deg,
y=1 * u.deg,
width=0.06 * u.deg,
length=0.3 * u.deg,
intensity=1500,
geometry=None,
):

if geometry is None:
s = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst")
geometry = s.tel[1].camera.geometry
geometry = s.tel[1].camera.geometry.transform_to(TelescopeFrame())

# make a toymodel shower model
model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi)
Expand Down Expand Up @@ -65,7 +62,7 @@ def test_hillas_selected(prod5_lst):
test Hillas-parameter routines on a sample image with selected values
against a sample image with masked values set to zero
"""
geom = prod5_lst.camera.geometry
geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame())
image, clean_mask = create_sample_image(geometry=geom)

image_zeros = image.copy()
Expand All @@ -81,15 +78,15 @@ def test_hillas_selected(prod5_lst):


def test_hillas_failure(prod5_lst):
geom = prod5_lst.camera.geometry
geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame())
blank_image = np.zeros(geom.n_pixels)

with pytest.raises(HillasParameterizationError):
hillas_parameters(geom, blank_image)


def test_hillas_masked_array(prod5_lst):
geom = prod5_lst.camera.geometry
geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame())
image, clean_mask = create_sample_image(psi="0d", geometry=geom)

image_zeros = image.copy()
Expand All @@ -103,30 +100,26 @@ def test_hillas_masked_array(prod5_lst):


def test_hillas_container(prod5_lst):
geom = prod5_lst.camera.geometry
geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame())
image, clean_mask = create_sample_image(psi="0d", geometry=geom)

params = hillas_parameters(geom[clean_mask], image[clean_mask])
assert isinstance(params, CameraHillasParametersContainer)

geom_telescope_frame = geom.transform_to(TelescopeFrame())
params = hillas_parameters(geom_telescope_frame[clean_mask], image[clean_mask])
assert isinstance(params, HillasParametersContainer)


def test_with_toy(prod5_lst):
rng = np.random.default_rng(42)

geom = prod5_lst.camera.geometry
geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame())

width = 0.03 * u.m
length = 0.15 * u.m
width_uncertainty = 0.00094 * u.m
length_uncertainty = 0.00465 * u.m
width = 0.06 * u.deg
length = 0.3 * u.deg
width_uncertainty = 0.0018 * u.deg
length_uncertainty = 0.0094 * u.deg
intensity = 500

xs = u.Quantity([0.5, 0.5, -0.5, -0.5], u.m)
ys = u.Quantity([0.5, -0.5, 0.5, -0.5], u.m)
xs = u.Quantity([1, 1, -1, -1], u.deg)
ys = u.Quantity([1, -1, 1, -1], u.deg)
psis = Angle([-90, -45, 0, 45, 90], unit="deg")

for x, y in zip(xs, ys):
Expand All @@ -141,8 +134,8 @@ def test_with_toy(prod5_lst):

result = hillas_parameters(geom, signal)

assert u.isclose(result.x, x, rtol=0.1)
assert u.isclose(result.y, y, rtol=0.1)
assert u.isclose(result.fov_lon, x, rtol=0.1)
assert u.isclose(result.fov_lat, y, rtol=0.1)

assert u.isclose(result.width, width, rtol=0.1)
assert u.isclose(result.width_uncertainty, width_uncertainty, rtol=0.4)
Expand All @@ -158,14 +151,14 @@ def test_with_toy(prod5_lst):
def test_skewness(prod5_lst):
rng = np.random.default_rng(42)

geom = prod5_lst.camera.geometry
geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame())

width = 0.03 * u.m
length = 0.15 * u.m
width = 0.10 * u.deg
length = 0.3 * u.deg
intensity = 2500

xs = u.Quantity([0.5, 0.5, -0.5, -0.5], u.m)
ys = u.Quantity([0.5, -0.5, 0.5, -0.5], u.m)
xs = u.Quantity([1, 1, -1, -1], u.deg)
ys = u.Quantity([1, -1, 1, -1], u.deg)
psis = Angle([-90, -45, 0, 45, 90], unit="deg")
skews = [0, 0.3, 0.6]

Expand All @@ -181,8 +174,8 @@ def test_skewness(prod5_lst):

result = hillas_parameters(geom, signal)

assert u.isclose(result.x, x, rtol=0.1)
assert u.isclose(result.y, y, rtol=0.1)
assert u.isclose(result.fov_lon, x, rtol=0.1)
assert u.isclose(result.fov_lat, y, rtol=0.1)

assert u.isclose(result.width, width, rtol=0.1)
assert u.isclose(result.length, length, rtol=0.1)
Expand All @@ -205,7 +198,7 @@ def test_skewness(prod5_lst):
def test_straight_line_width_0():
"""Test that hillas_parameters.width is 0 for a straight line of pixels"""
# three pixels in a straight line
long = np.array([0, 1, 2]) * 0.01
long = np.array([0, 1, 2]) * 1.25
trans = np.zeros(len(long))
pix_id = np.arange(len(long))

Expand All @@ -220,10 +213,11 @@ def test_straight_line_width_0():
geom = CameraGeometry(
name="testcam",
pix_id=pix_id,
pix_x=x * u.m,
pix_y=y * u.m,
pix_x=x * u.deg,
pix_y=y * u.deg,
pix_type="hexagonal",
pix_area=np.full(len(pix_id), 1.0) * u.m**2,
pix_area=np.full(len(pix_id), 1.0) * u.deg**2,
frame=TelescopeFrame(),
)

img = rng.poisson(5, size=len(long))
Expand All @@ -240,10 +234,11 @@ def test_single_pixel():
geom = CameraGeometry(
name="testcam",
pix_id=np.arange(9),
pix_x=x.ravel() * u.cm,
pix_y=y.ravel() * u.cm,
pix_x=x.ravel() * u.deg,
pix_y=y.ravel() * u.deg,
pix_type="rectangular",
pix_area=np.full(9, 1.0) * u.cm**2,
pix_area=np.full(9, 1.0) * u.deg**2,
frame=TelescopeFrame(),
)

image = np.zeros((3, 3))
Expand All @@ -255,82 +250,3 @@ def test_single_pixel():
assert hillas.length.value == 0
assert hillas.width.value == 0
assert np.isnan(hillas.psi)


def test_reconstruction_in_telescope_frame(prod5_lst):
"""
Compare the reconstruction in the telescope
and camera frame.
"""
np.random.seed(42)

geom = prod5_lst.camera.geometry
telescope_frame = TelescopeFrame()
camera_frame = geom.frame
geom_nom = geom.transform_to(telescope_frame)

width = 0.03 * u.m
length = 0.15 * u.m
intensity = 500

xs = u.Quantity([0.5, 0.5, -0.5, -0.5], u.m)
ys = u.Quantity([0.5, -0.5, 0.5, -0.5], u.m)
psis = Angle([-90, -45, 0, 45, 90], unit="deg")

def distance(coord):
return np.sqrt(np.diff(coord.x) ** 2 + np.diff(coord.y) ** 2) / 2

def get_transformed_length(telescope_hillas, telescope_frame, camera_frame):
main_edges = u.Quantity([-telescope_hillas.length, telescope_hillas.length])
main_lon = main_edges * np.cos(telescope_hillas.psi) + telescope_hillas.fov_lon
main_lat = main_edges * np.sin(telescope_hillas.psi) + telescope_hillas.fov_lat
cam_main_axis = SkyCoord(
fov_lon=main_lon, fov_lat=main_lat, frame=telescope_frame
).transform_to(camera_frame)
transformed_length = distance(cam_main_axis)
return transformed_length

def get_transformed_width(telescope_hillas, telescope_frame, camera_frame):
secondary_edges = u.Quantity([-telescope_hillas.width, telescope_hillas.width])
secondary_lon = (
secondary_edges * np.cos(telescope_hillas.psi) + telescope_result.fov_lon
)
secondary_lat = (
secondary_edges * np.sin(telescope_hillas.psi) + telescope_result.fov_lat
)
cam_secondary_edges = SkyCoord(
fov_lon=secondary_lon, fov_lat=secondary_lat, frame=telescope_frame
).transform_to(camera_frame)
transformed_width = distance(cam_secondary_edges)
return transformed_width

for x, y in zip(xs, ys):
for psi in psis:
# generate a toy image
model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi)
image, signal, noise = model.generate_image(
geom, intensity=intensity, nsb_level_pe=5
)

telescope_result = hillas_parameters(geom_nom, signal)
camera_result = hillas_parameters(geom, signal)
assert camera_result.intensity == telescope_result.intensity

# Compare results in both frames
transformed_cog = SkyCoord(
fov_lon=telescope_result.fov_lon,
fov_lat=telescope_result.fov_lat,
frame=telescope_frame,
).transform_to(camera_frame)
assert u.isclose(transformed_cog.x, camera_result.x, rtol=0.01)
assert u.isclose(transformed_cog.y, camera_result.y, rtol=0.01)

transformed_length = get_transformed_length(
telescope_result, telescope_frame, camera_frame
)
assert u.isclose(transformed_length, camera_result.length, rtol=0.01)

transformed_width = get_transformed_width(
telescope_result, telescope_frame, camera_frame
)
assert u.isclose(transformed_width, camera_result.width, rtol=0.01)
52 changes: 0 additions & 52 deletions ctapipe/image/tests/test_image_processor.py
Original file line number Diff line number Diff line change
@@ -1,17 +1,9 @@
"""
Tests for ImageProcessor functionality
"""
from copy import deepcopy

import astropy.units as u
import numpy as np
from numpy import isfinite

from ctapipe.calib import CameraCalibrator
from ctapipe.containers import (
CameraHillasParametersContainer,
CameraTimingParametersContainer,
)
from ctapipe.image import ImageProcessor
from ctapipe.image.cleaning import MARSImageCleaner

Expand Down Expand Up @@ -41,47 +33,3 @@ def test_image_processor(example_event, example_subarray):
assert isfinite(dl1tel.parameters.peak_time_statistics.max)

process_images.check_image.to_table()


def test_image_processor_camera_frame(example_event, example_subarray):
"""ensure we get parameters in the camera frame if explicitly specified"""
event = deepcopy(example_event)

calibrate = CameraCalibrator(subarray=example_subarray)
process_images = ImageProcessor(
subarray=example_subarray,
use_telescope_frame=False,
image_cleaner_type="MARSImageCleaner",
)

assert isinstance(process_images.clean, MARSImageCleaner)

calibrate(event)
process_images(event)

for dl1tel in event.dl1.tel.values():
assert isfinite(dl1tel.image_mask.sum())
assert isfinite(dl1tel.parameters.hillas.length.value)
dl1tel.parameters.hillas.length.to("meter")
assert isfinite(dl1tel.parameters.timing.slope.value)
assert isfinite(dl1tel.parameters.leakage.pixels_width_1)
assert isfinite(dl1tel.parameters.concentration.cog)
assert isfinite(dl1tel.parameters.morphology.n_pixels)
assert isfinite(dl1tel.parameters.intensity_statistics.max)
assert isfinite(dl1tel.parameters.peak_time_statistics.max)

process_images.check_image.to_table()

# set image to zeros to test invalid hillas parameters
# are in the correct frame
event = deepcopy(example_event)
calibrate(event)
for dl1 in event.dl1.tel.values():
dl1.image = np.zeros_like(dl1.image)

process_images(event)
for dl1 in event.dl1.tel.values():
assert isinstance(dl1.parameters.hillas, CameraHillasParametersContainer)
assert isinstance(dl1.parameters.timing, CameraTimingParametersContainer)
assert np.isnan(dl1.parameters.hillas.length.value)
assert dl1.parameters.hillas.length.unit == u.m
Loading

0 comments on commit 87073c6

Please sign in to comment.