Coverage for python / lsst / ip / isr / vignette.py: 23%

1# This file is part of ip_isr. 

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 

22__all__ = ('VignetteConfig', 'VignetteTask', 'setValidPolygonCcdIntersect', 

23 'maskVignettedRegion') 

24 

25import logging 

26import numpy as np 

27 

28import lsst.geom as geom 

29import lsst.afw.cameraGeom as cameraGeom 

30import lsst.afw.geom as afwGeom 

31 

32from lsst.pex.config import Config, Field 

33from lsst.pipe.base import Task 

34 

35 

36class VignetteConfig(Config): 

37 """Settings to define vignetting pattern. 

38 """ 

39 xCenter = Field( 

40 dtype=float, 

41 doc="Center of vignetting pattern, in focal plane x coordinates.", 

42 default=0.0, 

43 ) 

44 yCenter = Field( 

45 dtype=float, 

46 doc="Center of vignetting pattern, in focal plane y coordinates.", 

47 default=0.0, 

48 ) 

49 radius = Field( 

50 dtype=float, 

51 doc="Radius of vignetting pattern, in focal plane coordinates.", 

52 default=100.0, 

53 check=lambda x: x >= 0 

54 ) 

55 numPolygonPoints = Field( 

56 dtype=int, 

57 doc="Number of points used to define the vignette polygon.", 

58 default=100, 

59 ) 

60 

61 

62class VignetteTask(Task): 

63 """Define a simple circular vignette pattern and optionally update mask 

64 plane. 

65 """ 

66 ConfigClass = VignetteConfig 

67 _DefaultName = "isrVignette" 

68 

69 def run( 

70 self, 

71 exposure=None, 

72 doUpdateMask=True, 

73 maskPlane="NO_DATA", 

74 vignetteValue=None, 

75 log=None, 

76 doUpdatePolygon=True, 

77 ): 

78 """Generate circular vignette pattern. 

79 

80 Parameters 

81 ---------- 

82 exposure : `lsst.afw.image.Exposure`, optional 

83 Exposure to construct, apply, and optionally mask vignette for. 

84 doUpdateMask : `bool`, optional 

85 If true, the mask will be updated to mask the vignetted region. 

86 maskPlane : `str`, optional 

87 Mask plane to assign vignetted pixels to. 

88 vignetteValue : `float` or `None`, optional 

89 Value to assign to the image array pixels within the ``polygon`` 

90 region. If `None`, image pixel values are not replaced. 

91 log : `logging.Logger`, optional 

92 Log object to write to. 

93 doUpdatePolygon : `bool`, optional 

94 If true, the valid polygon will be updated. 

95 

96 Returns 

97 ------- 

98 polygon : `lsst.afw.geom.Polygon` 

99 Polygon defining the boundary of the vignetted region. 

100 """ 

101 theta = np.linspace(0, 2*np.pi, num=self.config.numPolygonPoints, endpoint=False) 

102 x = self.config.radius*np.cos(theta) + self.config.xCenter 

103 y = self.config.radius*np.sin(theta) + self.config.yCenter 

104 points = np.array([x, y]).transpose() 

105 fpPolygon = afwGeom.Polygon([geom.Point2D(x1, y1) for x1, y1 in reversed(points)]) 

106 if exposure is None: 

107 return fpPolygon 

108 

109 # Exposure was provided, so attach the validPolygon associated with the 

110 # vignetted region. 

111 validPolygon = setValidPolygonCcdIntersect(exposure, fpPolygon, log=log, setPolygon=doUpdatePolygon) 

112 

113 if doUpdateMask: 

114 maskVignettedRegion( 

115 exposure, 

116 validPolygon, 

117 maskPlane=maskPlane, 

118 vignetteValue=vignetteValue, 

119 log=log, 

120 ) 

121 

122 return fpPolygon 

123 

124 

125def setValidPolygonCcdIntersect(ccdExposure, fpPolygon, log=None, setPolygon=True): 

126 """Set valid polygon on ccdExposure associated with focal plane polygon. 

127 

128 The ccd exposure's valid polygon is the intersection of fpPolygon, 

129 a valid polygon in focal plane coordinates, and the ccd corners, 

130 in ccd pixel coordinates. 

131 

132 Parameters 

133 ---------- 

134 ccdExposure : `lsst.afw.image.Exposure` 

135 Exposure to process. 

136 fpPolygon : `lsst.afw.geom.Polygon` 

137 Polygon in focal plane coordinates. 

138 log : `logging.Logger`, optional 

139 Log object to write to. 

140 setPolygon : `bool`, optional 

141 Actually set the polygon, or just return it? 

142 

143 Returns 

144 ------- 

145 validPolygon : `lsst.afw.geom.Polygon` 

146 Valid polygon with the intersection. 

147 """ 

148 # Get ccd corners in focal plane coordinates 

149 ccd = ccdExposure.getDetector() 

150 fpCorners = ccd.getCorners(cameraGeom.FOCAL_PLANE) 

151 ccdPolygon = afwGeom.Polygon(fpCorners) 

152 # Get intersection of ccd corners with fpPolygon 

153 try: 

154 intersect = ccdPolygon.intersectionSingle(fpPolygon) 

155 except afwGeom.SinglePolygonException: 

156 intersect = None 

157 if intersect is not None: 

158 # Transform back to pixel positions and build new polygon 

159 ccdPoints = ccd.transform(intersect, cameraGeom.FOCAL_PLANE, cameraGeom.PIXELS) 

160 validPolygon = afwGeom.Polygon(ccdPoints) 

161 if setPolygon: 

162 ccdExposure.getInfo().setValidPolygon(validPolygon) 

163 else: 

164 if log is not None: 

165 log.info("Ccd exposure does not overlap with focal plane polygon. Not setting validPolygon.") 

166 validPolygon = None 

167 

168 return validPolygon 

169 

170 

171def maskVignettedRegion(exposure, polygon, maskPlane="NO_DATA", vignetteValue=None, log=None): 

172 """Add mask bit to image pixels according to vignetted polygon region. 

173 

174 NOTE: this function could be used to mask and replace pixels associated 

175 with any polygon in the exposure pixel coordinates. 

176 

177 Parameters 

178 ---------- 

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

180 Image whose mask plane is to be updated. 

181 polygon : `lsst.afw.geom.Polygon` 

182 Polygon region defining the vignetted region in the pixel coordinates 

183 of ``exposure``. 

184 maskPlane : `str`, optional 

185 Mask plane to assign vignetted pixels to. 

186 vignetteValue : `float` or `None`, optional 

187 Value to assign to the image array pixels within the ``polygon`` 

188 region. If `None`, image pixel values are not replaced. 

189 log : `logging.Logger`, optional 

190 Log object to write to. 

191 

192 Raises 

193 ------ 

194 RuntimeError 

195 Raised if no valid polygon exists. 

196 """ 

197 log = log if log else logging.getLogger(__name__) 

198 if not polygon: 

199 log.info("No polygon provided. Masking nothing.") 

200 return 

201 

202 fullyIlluminated = True 

203 if not all(polygon.contains(exposure.getBBox().getCorners())): 

204 fullyIlluminated = False 

205 log.info("Exposure is fully illuminated? %s", fullyIlluminated) 

206 

207 if not fullyIlluminated: 

208 # Scan pixels. 

209 mask = exposure.getMask() 

210 xx, yy = np.meshgrid(np.arange(0, mask.getWidth(), dtype=float), 

211 np.arange(0, mask.getHeight(), dtype=float)) 

212 

213 vignMask = ~(polygon.contains(xx, yy)) 

214 

215 bitMask = mask.getPlaneBitMask(maskPlane) 

216 maskArray = mask.getArray() 

217 maskArray[vignMask] |= bitMask 

218 log.info("Exposure contains {} vignetted pixels which are now masked with mask plane {}.". 

219 format(np.count_nonzero(vignMask), maskPlane)) 

220 

221 if vignetteValue is not None: 

222 imageArray = exposure.getImage().getArray() 

223 imageArray[vignMask] = vignetteValue 

224 log.info("Vignetted pixels in image array have been replaced with {}.".format(vignetteValue))