Coverage for rta_reconstruction/image_displacement.py: 28%

1import astropy.units as u 

2import numpy as np 

3 

4from rta_reconstruction.utils.coordinates import pol2cart, sky_to_camera 

5 

6 

7def disp(cog_x, cog_y, src_x, src_y, hillas_psi): 

8 """ 

9 Compute the disp parameters 

10 

11 Parameters 

12 ---------- 

13 cog_x: `numpy.ndarray` or float 

14 cog_y: `numpy.ndarray` or float 

15 src_x: `numpy.ndarray` or float 

16 src_y: `numpy.ndarray` or float 

17 hillas_psi: `numpy.ndarray` or float 

18 

19 Returns 

20 ------- 

21 (disp_dx, disp_dy, disp_norm, disp_angle, disp_sign): 

22 disp_dx: 'astropy.units.m` 

23 disp_dy: 'astropy.units.m` 

24 disp_norm: 'astropy.units.m` 

25 disp_angle: 'astropy.units.rad` 

26 disp_sign: `numpy.ndarray` 

27 """ 

28 disp_dx = src_x - cog_x 

29 disp_dy = src_y - cog_y 

30 

31 disp_norm = disp_dx * np.cos(hillas_psi) + disp_dy * np.sin(hillas_psi) 

32 disp_sign = np.sign(disp_norm) 

33 disp_norm = np.abs(disp_norm) 

34 

35 # disp_sign : indicates in which direction, "positive" or "negative", we must move along the 

36 # reconstructed image axis (with direction defined by the versor cos(hillas_psi), sin(hillas_psi)) 

37 # we must move from cog_x, cog_y to get closest to the true direction (src_x, src_y) 

38 

39 if hasattr(disp_dx, "__len__"): 

40 disp_angle = np.arctan(disp_dy / disp_dx) 

41 disp_angle[disp_dx == 0] = np.pi / 2.0 * np.sign(disp_dy[disp_dx == 0]) 

42 else: 

43 if disp_dx == 0: 

44 disp_angle = np.pi / 2.0 * np.sign(disp_dy) 

45 else: 

46 disp_angle = np.arctan(disp_dy / disp_dx) 

47 

48 return disp_dx, disp_dy, disp_norm, disp_angle, disp_sign 

49 

50 

51def miss(disp_dx, disp_dy, hillas_psi): 

52 """ 

53 Compute miss 

54 

55 Parameters 

56 ---------- 

57 disp_dx: `numpy.ndarray` or float 

58 disp_dy: `numpy.ndarray` or float 

59 hillas_psi: `numpy.ndarray` or float 

60 

61 Returns 

62 ------- 

63 

64 """ 

65 return np.abs(np.sin(hillas_psi) * disp_dx - np.cos(hillas_psi) * disp_dy) 

66 

67 

68def disp_vector_pol2cart(disp_norm, disp_angle, disp_sign): 

69 """ 

70 Compute `disp_norm.dx` and `disp_norm.dy` vector from `disp_norm.norm`, 

71 `disp_norm.angle` and `disp_norm.sign` 

72 

73 Parameters 

74 ---------- 

75 disp_norm: float 

76 disp_angle: float 

77 disp_sign: float 

78 

79 Returns 

80 ------- 

81 disp_dx, disp_dy 

82 """ 

83 return np.transpose(pol2cart(disp_norm, disp_angle, disp_sign)) 

84 

85 

86def disp_to_pos(disp_dx, disp_dy, cog_x, cog_y): 

87 """ 

88 Calculates source position in camera coordinates(x,y) 

89 from the reconstructed disp. 

90 

91 Parameters 

92 ---------- 

93 disp: DispContainer 

94 cog_x: float 

95 Coordinate x of the center of gravity of Hillas ellipse 

96 cog_y: float 

97 Coordinate y of the center of gravity of Hillas ellipse 

98 

99 Returns 

100 ------- 

101 (source_pos_x, source_pos_y) 

102 """ 

103 source_pos_x = cog_x + disp_dx 

104 source_pos_y = cog_y + disp_dy 

105 

106 return source_pos_x, source_pos_y 

107 

108 

109def update_disp_with_effective_focal_length(data, effective_focal_length=29.30565 * u.m): 

110 """Update disp parameters using effective focal length 

111 

112 Parameters 

113 ---------- 

114 data: Pandas DataFrame 

115 config: dictionnary containing configuration 

116 """ 

117 

118 source_pos_in_camera = sky_to_camera( 

119 u.Quantity(data["mc_alt"].values, u.rad, copy=False), 

120 u.Quantity(data["mc_az"].values, u.rad, copy=False), 

121 effective_focal_length, 

122 u.Quantity(data["mc_alt_tel"].values, u.rad, copy=False), 

123 u.Quantity(data["mc_az_tel"].values, u.rad, copy=False), 

124 ) 

125 

126 expected_src_pos_x_m = source_pos_in_camera.x.to_value(u.m) 

127 expected_src_pos_y_m = source_pos_in_camera.y.to_value(u.m) 

128 

129 data["src_x"] = expected_src_pos_x_m 

130 data["src_y"] = expected_src_pos_y_m 

131 

132 disp_dx, disp_dy, disp_norm, disp_angle, disp_sign = disp.disp( 

133 data["x"].values, data["y"].values, expected_src_pos_x_m, expected_src_pos_y_m, data["psi"].values 

134 ) 

135 

136 data["disp_dx"] = disp_dx 

137 data["disp_dy"] = disp_dy 

138 data["disp_norm"] = disp_norm 

139 data["disp_angle"] = disp_angle 

140 data["disp_sign"] = disp_sign