diff --git a/ctapipe/image/extractor.py b/ctapipe/image/extractor.py index 2cf3eef5519..91215f5d132 100644 --- a/ctapipe/image/extractor.py +++ b/ctapipe/image/extractor.py @@ -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 diff --git a/ctapipe/image/tests/test_concentration.py b/ctapipe/image/tests/test_concentration.py index 25dd14e0501..c0ed4241fb8 100644 --- a/ctapipe/image/tests/test_concentration.py +++ b/ctapipe/image/tests/test_concentration.py @@ -1,6 +1,7 @@ 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 @@ -8,7 +9,7 @@ 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]) @@ -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 diff --git a/ctapipe/image/tests/test_hillas.py b/ctapipe/image/tests/test_hillas.py index b3b148c8bc6..35e981babaa 100644 --- a/ctapipe/image/tests/test_hillas.py +++ b/ctapipe/image/tests/test_hillas.py @@ -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 @@ -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) @@ -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() @@ -81,7 +78,7 @@ 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): @@ -89,7 +86,7 @@ def test_hillas_failure(prod5_lst): 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() @@ -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): @@ -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) @@ -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] @@ -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) @@ -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)) @@ -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)) @@ -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)) @@ -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) diff --git a/ctapipe/image/tests/test_image_processor.py b/ctapipe/image/tests/test_image_processor.py index 327b17aebf6..63c50265940 100644 --- a/ctapipe/image/tests/test_image_processor.py +++ b/ctapipe/image/tests/test_image_processor.py @@ -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 @@ -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 diff --git a/ctapipe/image/tests/test_timing_parameters.py b/ctapipe/image/tests/test_timing_parameters.py index a37c07b12db..9dd3bb9abb4 100644 --- a/ctapipe/image/tests/test_timing_parameters.py +++ b/ctapipe/image/tests/test_timing_parameters.py @@ -2,7 +2,8 @@ import numpy as np from numpy.testing import assert_allclose -from ctapipe.containers import CameraHillasParametersContainer +from ctapipe.containers import HillasParametersContainer +from ctapipe.coordinates import TelescopeFrame def test_psi_0(prod5_lst): @@ -16,8 +17,10 @@ def test_psi_0(prod5_lst): intercept = 1.0 deviation = 0.1 - geom = prod5_lst.camera.geometry - hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=0 * u.deg) + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) + hillas = HillasParametersContainer( + fov_lon=0 * u.deg, fov_lat=0 * u.deg, psi=0 * u.deg + ) random = np.random.default_rng(0) peak_time = intercept + grad * geom.pix_x.value @@ -45,9 +48,9 @@ def test_psi_20(prod5_lst): intercept = 1 deviation = 0.1 - geom = prod5_lst.camera.geometry + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) psi = 20 * u.deg - hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=psi) + hillas = HillasParametersContainer(fov_lon=0 * u.deg, fov_lat=0 * u.deg, psi=psi) random = np.random.default_rng(0) peak_time = intercept + grad * ( @@ -76,8 +79,10 @@ def test_ignore_negative(prod5_lst): intercept = 1.0 deviation = 0.1 - geom = prod5_lst.camera.geometry - hillas = CameraHillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=0 * u.deg) + geom = prod5_lst.camera.geometry.transform_to(TelescopeFrame()) + hillas = HillasParametersContainer( + fov_lon=0 * u.deg, fov_lat=0 * u.deg, psi=0 * u.deg + ) random = np.random.default_rng(0) peak_time = intercept + grad * geom.pix_x.value