Coverage for tests/test_findGlintTrails.py: 93%

1# This file is part of meas_algorithms. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22import unittest 

23 

24import numpy as np 

25 

26from lsst.meas.algorithms import findGlintTrails 

27import lsst.geom 

28import lsst.meas.base.tests 

29import lsst.utils.tests 

30 

31# Set this and have display_ds9 setup to see the fit on the image. 

32display = False 

33if display: 33 ↛ 34line 33 didn't jump to line 34 because the condition on line 33 was never true

34 import lsst.afw.display 

35 display = lsst.afw.display.Display() 

36 

37 

38class TestFindGlintTrails(lsst.utils.tests.TestCase): 

39 """Generate two glint trails and test that they are both found, while 

40 outlier points are not included. 

41 """ 

42 

43 def _make_image_3_trails(self): 

44 """Make an image with 3 trails on it, one of them negative. 

45 """ 

46 rng = np.random.default_rng(10) 

47 bbox = lsst.geom.Box2I(lsst.geom.Point2I(-5, -4), lsst.geom.Point2I(1005, 1084)) 

48 dataset = lsst.meas.base.tests.TestDataset(bbox) 

49 flux = 10000 

50 

51 x0 = 200 

52 y0 = 300 

53 scale = 50 

54 # one outlier source near the first glint trail. 

55 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(x0 - 31, y0 - 27)) 

56 # two glint trails 

57 for i in range(8): 

58 dataset.addSource(instFlux=flux, 

59 centroid=lsst.geom.Point2D(i*scale + x0 + rng.random(), 

60 i*scale + y0 + rng.random())) 

61 

62 # A negative trail (e.g. came from the template) that shouldn't be 

63 # included in any of the fits. 

64 x0 = 170 

65 y0 = 690 

66 step = 70 

67 for i in range(6): 

68 dataset.addSource(instFlux=-flux, 

69 centroid=lsst.geom.Point2D(i*step + x0 + rng.random(), 

70 -i*step + y0 + rng.random()), 

71 negative=True) 

72 

73 # a source that shouldn't appear in the found trails 

74 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(500, 200)) 

75 

76 # Perpendicular trail that has a point lying on the same line as the 

77 # first trail. 

78 x0 = 380 

79 y0 = 760 

80 step = 70 

81 for i in range(5): 

82 dataset.addSource(instFlux=flux, 

83 centroid=lsst.geom.Point2D(i*step + x0 + rng.random(), 

84 -i*step + y0 + rng.random())) 

85 

86 schema = dataset.makeMinimalSchema() 

87 schema.addField("ip_diffim_DipoleFit_classification", type="Flag") 

88 exposure, catalog = dataset.realize(10.0, schema=schema) 

89 # So that the catalog ids don't look like simple indices when debugging. 

90 catalog["id"] += 100 

91 return exposure, catalog 

92 

93 def test_simple(self): 

94 """Test the basic operation of the glint finder on an image with two 

95 glint trails. 

96 """ 

97 exposure, catalog = self._make_image_3_trails() 

98 

99 config = findGlintTrails.FindGlintTrailsTask.ConfigClass() 

100 # Limit the search radius so that we can test the line-extension 

101 # part of the fitter. 

102 config.radius = 300 

103 # Use a tighter threshold for this simulated data, which has better 

104 # controlled centroids than real data. 

105 config.threshold = 1.5 

106 task = findGlintTrails.FindGlintTrailsTask(config=config) 

107 result = task.run(catalog) 

108 

109 if display: 109 ↛ 110line 109 didn't jump to line 110 because the condition on line 109 was never true

110 display.frame = 1 

111 display.image(exposure, title="something") 

112 for i, trail in enumerate(result.trails): 

113 display.centroids(trail, size=5, ctype="cyan", symbol=i) 

114 

115 self.assertEqual(len(result.trails), 2) 

116 # Note that if you add sources to the simulated catalog, the ids in 

117 # the fitted trails may change. 

118 self.assertSetEqual(set(result.trails[0]["id"]), 

119 {102, 103, 104, 105, 106, 107, 108, 109, 119}) 

120 self.assertSetEqual(set(result.trails[1]["id"]), 

121 {117, 118, 119, 120, 121}) 

122 self.assertSetEqual(result.trailed_ids, {102, 103, 104, 105, 106, 107, 108, 

123 109, 117, 118, 119, 120, 121}) 

124 # Expected lengths from the step size of the simulated trails. 

125 self.assertFloatsAlmostEqual(result.parameters[0].length, 

126 np.sqrt((7*50)**2 + (7*50)**2), 

127 rtol=1e-3) 

128 self.assertFloatsAlmostEqual(result.parameters[1].length, 

129 np.sqrt((4*70)**2 + (4*70)**2), 

130 rtol=1e-3) 

131 self.assertFloatsAlmostEqual(np.degrees(result.parameters[0].angle), 45.0, rtol=2e-3) 

132 self.assertFloatsAlmostEqual(np.degrees(result.parameters[1].angle), -45.0, rtol=2e-3) 

133 

134 def _make_catalog(self, centroids): 

135 """Make a catalog with one source at each of the given centroids.""" 

136 bbox = lsst.geom.Box2I(lsst.geom.Point2I(-5, -4), lsst.geom.Point2I(1005, 1084)) 

137 dataset = lsst.meas.base.tests.TestDataset(bbox) 

138 for x, y in centroids: 

139 dataset.addSource(instFlux=10000, centroid=lsst.geom.Point2D(x, y)) 

140 schema = dataset.makeMinimalSchema() 

141 schema.addField("ip_diffim_DipoleFit_classification", type="Flag") 

142 exposure, catalog = dataset.realize(10.0, schema=schema) 

143 return catalog 

144 

145 def test_trail_angles(self): 

146 """Test that trails at any angle, including vertical, are found and 

147 measured, and that a source further than ``threshold`` from the line 

148 is left out. 

149 """ 

150 rng = np.random.default_rng(20) 

151 step = 60 

152 task = findGlintTrails.FindGlintTrailsTask() 

153 for angle in (0, 30, 60, 80, 90, 135): 

154 with self.subTest(angle=angle): 

155 direction = np.array([np.cos(np.radians(angle)), np.sin(np.radians(angle))]) 

156 normal = np.array([-direction[1], direction[0]]) 

157 start = np.array([300.0, 300.0]) 

158 trail = [start + i*step*direction + 0.5*rng.random(2) for i in range(6)] 

159 # Within the 2x threshold candidate search, beyond threshold. 

160 outlier = start + 2.5*step*direction + 1.5*task.config.threshold*normal 

161 catalog = self._make_catalog(trail + [outlier]) 

162 

163 result = task.run(catalog) 

164 

165 self.assertEqual(len(result.trails), 1) 

166 # Source ids count from 1 in the order the sources were added. 

167 self.assertSetEqual(set(result.trails[0]["id"]), set(range(1, 7))) 

168 parameters = result.parameters[0] 

169 # Compare angles modulo 180 degrees: a vertical line is equally 

170 # well described by +90 and -90 degrees. 

171 offset = (np.degrees(parameters.angle) - angle + 90) % 180 - 90 

172 self.assertLess(abs(offset), 0.5) 

173 self.assertFloatsAlmostEqual(parameters.length, 5*step, atol=1.0) 

174 center = np.mean(trail, axis=0) 

175 self.assertFloatsAlmostEqual(parameters.center_x, center[0], atol=1e-6) 

176 self.assertFloatsAlmostEqual(parameters.center_y, center[1], atol=1e-6) 

177 fit_normal = np.array([-np.sin(parameters.angle), np.cos(parameters.angle)]) 

178 distances = (np.array(trail) - center) @ fit_normal 

179 self.assertFloatsAlmostEqual(parameters.stderr, np.sqrt(np.mean(distances**2)), 

180 rtol=1e-8) 

181 

182 def test_trail_of_min_points(self): 

183 """Test that a trail of exactly ``min_points`` sources, with no other 

184 sources nearby, is found. 

185 """ 

186 task = findGlintTrails.FindGlintTrailsTask() 

187 direction = np.array([np.cos(np.radians(30)), np.sin(np.radians(30))]) 

188 trail = [np.array([300.0, 300.0]) + i*60*direction for i in range(task.config.min_points)] 

189 catalog = self._make_catalog(trail) 

190 

191 result = task.run(catalog) 

192 

193 self.assertEqual(len(result.trails), 1) 

194 self.assertEqual(len(result.trails[0]), task.config.min_points) 

195 

196 def test_fit_line(self): 

197 """Test the perpendicular-distance line fit on a vertical line.""" 

198 x = np.array([10.0, 11.0, 10.0, 9.0]) 

199 y = np.array([0.0, 10.0, 20.0, 30.0]) 

200 

201 line = findGlintTrails.fit_line(x, y) 

202 

203 self.assertFloatsAlmostEqual(line.center, np.array([10.0, 15.0]), atol=1e-12) 

204 self.assertFloatsAlmostEqual(abs(line.direction[1]), 1.0, atol=1e-2) 

205 self.assertFloatsAlmostEqual(line.distance(13.0, 15.0), 3.0, atol=0.1) 

206 

207 def test_empty_image(self): 

208 """Test that no trails are found on an empty image.""" 

209 bbox = lsst.geom.Box2I(lsst.geom.Point2I(-5, -4), lsst.geom.Point2I(1005, 1084)) 

210 dataset = lsst.meas.base.tests.TestDataset(bbox) 

211 schema = dataset.makeMinimalSchema() 

212 schema.addField("ip_diffim_DipoleFit_classification", type="Flag") 

213 exposure, catalog = dataset.realize(10.0, schema=schema) 

214 

215 task = findGlintTrails.FindGlintTrailsTask() 

216 result = task.run(catalog) 

217 

218 self.assertEqual(len(result.trails), 0) 

219 self.assertEqual(result.trailed_ids, set()) 

220 

221 def test_many_points_no_trail(self): 

222 """Test that no trails are found on an image with many points, but no 

223 sets of three that could lie on the same line. 

224 """ 

225 bbox = lsst.geom.Box2I(lsst.geom.Point2I(-5, -4), lsst.geom.Point2I(1005, 1084)) 

226 dataset = lsst.meas.base.tests.TestDataset(bbox) 

227 flux = 10000 

228 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(500, 500)) 

229 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(480, 610)) 

230 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(520, 680)) 

231 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(570, 500)) 

232 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(640, 630)) 

233 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(600, 700)) 

234 schema = dataset.makeMinimalSchema() 

235 schema.addField("ip_diffim_DipoleFit_classification", type="Flag") 

236 exposure, catalog = dataset.realize(10.0, schema=schema) 

237 

238 task = findGlintTrails.FindGlintTrailsTask() 

239 result = task.run(catalog) 

240 

241 self.assertEqual(len(result.trails), 0) 

242 self.assertEqual(result.trailed_ids, set()) 

243 

244 def test_rejected_trail_min_points(self): 

245 """Test that no trails are found on an image with an initially 

246 reasonable trail that is rejected during fitting for having too many 

247 outliers. 

248 """ 

249 bbox = lsst.geom.Box2I(lsst.geom.Point2I(-5, -4), lsst.geom.Point2I(1005, 1084)) 

250 dataset = lsst.meas.base.tests.TestDataset(bbox) 

251 flux = 10000 

252 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(100, 100)) 

253 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(210, 200)) 

254 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(300, 310)) 

255 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(400, 400)) 

256 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(500, 500)) 

257 dataset.addSource(instFlux=flux, centroid=lsst.geom.Point2D(600, 600)) 

258 schema = dataset.makeMinimalSchema() 

259 schema.addField("ip_diffim_DipoleFit_classification", type="Flag") 

260 exposure, catalog = dataset.realize(10.0, schema=schema) 

261 

262 task = findGlintTrails.FindGlintTrailsTask() 

263 result = task.run(catalog) 

264 

265 self.assertEqual(len(result.trails), 0) 

266 self.assertEqual(result.trailed_ids, set()) 

267 

268 

269def setup_module(module): 

270 lsst.utils.tests.init() 

271 

272 

273if __name__ == "__main__": 273 ↛ 274line 273 didn't jump to line 274 because the condition on line 273 was never true

274 lsst.utils.tests.init() 

275 unittest.main()