Coverage for python / lsst / obs / subaru / strayLight / yStrayLight.py: 20%

1# Copyright (C) 2017 HSC Software Team 

2# Copyright (C) 2017 Sogo Mineo 

3# 

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

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

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

7# (at your option) any later version. 

8# 

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

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

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

12# GNU General Public License for more details. 

13# 

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

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

16 

17__all__ = ["SubaruStrayLightTask"] 

18 

19from typing import Optional 

20 

21import numpy 

22from astropy.io import fits 

23import scipy.interpolate 

24 

25from lsst.daf.base import DateTime 

26from lsst.geom import Angle, degrees 

27from lsst.daf.butler import DeferredDatasetHandle 

28from lsst.ip.isr.straylight import StrayLightConfig, StrayLightTask, StrayLightData 

29 

30from . import waveletCompression 

31from .rotatorAngle import inrStartEnd 

32 

33 

34BAD_THRESHOLD = 500 # Threshold for identifying bad pixels in the reconstructed dark image 

35 

36# LEDs causing the stray light were been covered up by his date. 

37# We believe there is no remaining stray light. 

38_STRAY_LIGHT_FIXED_DATE = DateTime("2018-01-01T00:00:00Z", DateTime.UTC) 

39 

40 

41class SubaruStrayLightTask(StrayLightTask): 

42 """Remove stray light in the y-band 

43 

44 LEDs in an encoder in HSC are producing stray light on the detectors, 

45 producing the "Eye of Y-band" feature. It can be removed by 

46 subtracting open-shutter darks. However, because the pattern of stray 

47 light varies with rotator angle, many dark exposures are required. 

48 To reduce the data volume for the darks, the images have been 

49 compressed using wavelets. The code used to construct these is at: 

50 

51 https://hsc-gitlab.mtk.nao.ac.jp/sogo.mineo/ybackground/ 

52 

53 This Task retrieves the appropriate dark, uncompresses it and 

54 uses it to remove the stray light from an exposure. 

55 """ 

56 

57 ConfigClass = StrayLightConfig 

58 

59 def check(self, exposure): 

60 # Docstring inherited from StrayLightTask.check. 

61 detId = exposure.getDetector().getId() 

62 if not self.checkFilter(exposure): 

63 # No correction to be made 

64 return False 

65 if detId in range(104, 112): 

66 # No correction data: assume it's zero 

67 return False 

68 if exposure.visitInfo.date >= _STRAY_LIGHT_FIXED_DATE: 

69 return False 

70 

71 return True 

72 

73 def run(self, exposure, strayLightData): 

74 """Subtract the y-band stray light 

75 

76 This relies on knowing the instrument rotator angle during the 

77 exposure. The FITS headers provide only the instrument rotator 

78 angle at the start of the exposure (INR_STR), but better 

79 stray light removal is obtained when we calculate the start and 

80 end instrument rotator angles ourselves (config parameter 

81 ``doRotatorAngleCorrection=True``). 

82 

83 Parameters 

84 ---------- 

85 exposure : `lsst.afw.image.Exposure` 

86 Exposure to correct. 

87 strayLightData : `SubaruStrayLightData` or 

88 `~lsst.daf.butler.DeferredDatasetHandle` 

89 An opaque object that contains any calibration data used to 

90 correct for stray light. 

91 """ 

92 if not self.check(exposure): 

93 return None 

94 

95 if strayLightData is None: 

96 raise RuntimeError("No strayLightData supplied for correction.") 

97 

98 if isinstance(strayLightData, DeferredDatasetHandle): 

99 # Get the deferred object. 

100 strayLightData = strayLightData.get() 

101 

102 exposureMetadata = exposure.getMetadata() 

103 detId = exposure.getDetector().getId() 

104 if self.config.doRotatorAngleCorrection: 

105 angleStart, angleEnd = inrStartEnd(exposure.getInfo().getVisitInfo()) 

106 self.log.debug( 

107 "(INR-STR, INR-END) = (%g, %g) (FITS header says (%g, %g)).", 

108 angleStart, angleEnd, 

109 exposureMetadata.getDouble('INR-STR'), exposureMetadata.getDouble('INR-END') 

110 ) 

111 else: 

112 angleStart = exposureMetadata.getDouble('INR-STR') 

113 angleEnd = None 

114 

115 self.log.info("Correcting y-band background.") 

116 

117 model = strayLightData.evaluate(angleStart*degrees, 

118 None if angleStart == angleEnd else angleEnd*degrees) 

119 

120 # Some regions don't have useful model values because the amplifier is 

121 # dead when the darks were taken 

122 # 

123 badAmps = {9: [0, 1, 2, 3], 33: [0, 1], 43: [0]} # Known bad amplifiers in the data: {ccdId: [ampId]} 

124 if detId in badAmps: 

125 isBad = numpy.zeros_like(model, dtype=bool) 

126 for ii in badAmps[detId]: 

127 amp = exposure.getDetector()[ii] 

128 box = amp.getBBox() 

129 isBad[box.getBeginY():box.getEndY(), box.getBeginX():box.getEndX()] = True 

130 mask = exposure.getMaskedImage().getMask() 

131 if numpy.all(isBad): 

132 model[:] = 0.0 

133 else: 

134 model[isBad] = numpy.median(model[~isBad]) 

135 mask.array[isBad] |= mask.getPlaneBitMask("SUSPECT") 

136 

137 model *= exposure.getInfo().getVisitInfo().getExposureTime() 

138 exposure.image.array -= model 

139 

140 

141class SubaruStrayLightData(StrayLightData): 

142 """Object that reads and integrates the wavelet-compressed 

143 HSC y-band stray-light model. 

144 

145 Parameters 

146 ---------- 

147 filename : `str` 

148 Full path to a FITS files containing the stray-light model. 

149 """ 

150 

151 @classmethod 

152 def readFits(cls, filename, **kwargs): 

153 calib = cls() 

154 

155 with fits.open(filename) as hdulist: 

156 calib.ampData = [hdu.data for hdu in hdulist] 

157 calib.setMetadata(hdulist[0].header) 

158 

159 calib.log.info("Finished reading straylightData.") 

160 return calib 

161 

162 def evaluate(self, angle_start: Angle, angle_end: Optional[Angle] = None): 

163 """Get y-band background image array for a range of angles. 

164 

165 It is hypothesized that the instrument rotator rotates at a constant 

166 angular velocity. This is not strictly true, but should be a 

167 sufficient approximation for the relatively short exposure times 

168 typical for HSC. 

169 

170 Parameters 

171 ---------- 

172 angle_start : `float` 

173 Instrument rotation angle in degrees at the start of the exposure. 

174 angle_end : `float`, optional 

175 Instrument rotation angle in degrees at the end of the exposure. 

176 If not provided, the returned array will reflect a snapshot at 

177 `angle_start`. 

178 

179 Returns 

180 ------- 

181 ccd_img : `numpy.ndarray` 

182 Background data for this exposure. 

183 """ 

184 header = self.getMetadata() 

185 

186 # full-size ccd height & channel width 

187 ccd_h, ch_w = header["F_NAXIS2"], header["F_NAXIS1"] 

188 # saved data is compressed to 1/2**scale_level of the original size 

189 image_scale_level = header["WTLEVEL2"], header["WTLEVEL1"] 

190 angle_scale_level = header["WTLEVEL3"] 

191 

192 ccd_w = ch_w * len(self.ampData) 

193 ccd_img = numpy.empty(shape=(ccd_h, ccd_w), dtype=numpy.float32) 

194 

195 for ch, hdu in enumerate(self.ampData): 

196 volume = _upscale_volume(hdu, angle_scale_level) 

197 

198 if angle_end is None: 

199 img = volume(angle_start.asDegrees()) 

200 else: 

201 img = (volume.integrate(angle_start.asDegrees(), angle_end.asDegrees()) 

202 * (1.0 / (angle_end.asDegrees() - angle_start.asDegrees()))) 

203 

204 ccd_img[:, ch_w*ch:ch_w*(ch+1)] = _upscale_image(img, (ccd_h, ch_w), image_scale_level) 

205 

206 # Some regions don't have useful values because the amplifier is dead 

207 # when the darks were taken 

208 # is_bad = ccd_img > BAD_THRESHOLD 

209 # ccd_img[is_bad] = numpy.median(ccd_img[~is_bad]) 

210 

211 return ccd_img 

212 

213 

214def _upscale_image(img, target_shape, level): 

215 """Upscale the given image to `target_shape` . 

216 

217 Parameters 

218 ---------- 

219 img : `numpy.array`, (Nx, Ny) 

220 Compressed image. ``img.shape`` must agree 

221 with waveletCompression.scaled_size(target_shape, level) 

222 target_shape : `tuple` [`int`, `int`] 

223 The shape of upscaled image, which is to be returned. 

224 level : `int` or `tuple` [`int`] 

225 Level of multiresolution analysis (or synthesis) 

226 

227 Returns 

228 ------- 

229 resized : `numpy.array`, (Nu, Nv) 

230 Upscaled image with the ``target_shape``. 

231 """ 

232 h, w = waveletCompression.scaled_size(target_shape, level) 

233 

234 large_img = numpy.zeros(shape=target_shape, dtype=float) 

235 large_img[:h, :w] = img 

236 

237 return waveletCompression.icdf_9_7(large_img, level) 

238 

239 

240def _upscale_volume(volume, level): 

241 """Upscale the given volume (= sequence of images) along the 0-th 

242 axis, and return an instance of a interpolation object that 

243 interpolates the 0-th axis. The 0-th axis is the instrument 

244 rotation. 

245 

246 Parameters 

247 ---------- 

248 volume : `numpy.array`, (Nx, Ny, Nz) 

249 Sequence of images. 

250 level : `int` 

251 Level of multiresolution analysis along the 0-th axis. 

252 

253 interpolator : callable 

254 An object that returns a slice of the volume at a specific 

255 angle (in degrees), with one positional argument: 

256 

257 - ``angle``: The angle in degrees. 

258 """ 

259 angles = 720 

260 _, h, w = volume.shape 

261 

262 large_volume = numpy.zeros(shape=(angles+1, h, w), dtype=float) 

263 

264 layers = waveletCompression.scaled_size(angles, level) 

265 large_volume[:layers] = volume 

266 

267 large_volume[:-1] = waveletCompression.periodic_icdf_9_7_1d(large_volume[:-1], level, axis=0) 

268 large_volume[-1] = large_volume[0] 

269 

270 x = numpy.arange(angles+1) / 2.0 

271 return scipy.interpolate.CubicSpline(x, large_volume, axis=0, bc_type="periodic")