Coverage for python / lsst / obs / subaru / crosstalkYagi.py: 0%
108 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-22 09:23 +0000
1#
2# LSST Data Management System
3# Copyright 2008, 2009, 2010, 2011 LSST Corporation.
4#
5# This product includes software developed by the
6# LSST Project (http://www.lsst.org/).
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <http://www.lsstcorp.org/LegalNotices/>.
21#
22"""
23Determine and apply crosstalk corrections
25N.b. This code was written and tested for the 4-amplifier Hamamatsu chips used
26in (Hyper)?SuprimeCam, and will need to be generalised to handle other
27amplifier layouts. I don't want to do this until we have an example.
29N.b. To estimate crosstalk from the SuprimeCam data, the commands are e.g.:
31.. code-block:: python
33 import crosstalk
34 coeffs = crosstalk.estimateCoeffs(range(131634, 131642), range(10),
35 threshold=1e5, plot=True, title="CCD0..9",
36 fig=1)
37 crosstalk.fixCcd(131634, 0, coeffs)
38"""
39import numpy as np
40import lsst.afw.detection as afwDetect
41import lsst.afw.image as afwImage
42import lsst.afw.math as afwMath
43import lsst.afw.display as afwDisplay
44import lsst.geom as geom
45import lsst.pipe.base as pipeBase
46import lsst.pex.config as pexConfig
47from lsst.obs.subaru import subtractCrosstalk
48from lsst.obs.subaru.crosstalk import makeList, estimateCoeffs, printCoeffs, readImage, subtractXTalk
51class CrosstalkYagiCoeffsConfig(pexConfig.Config):
52 """Specify crosstalk coefficients for a CCD"""
54 crossTalkCoeffs1 = pexConfig.ListField(
55 dtype=float,
56 doc="crosstalk coefficients (primary by the high-count pixel)",
57 default=[0, -0.000148, -0.000162, -0.000167, # cAA,cAB,cAC,cAD
58 -0.000148, 0, -0.000077, -0.000162, # cBA,cBB,cBC,cBD
59 -0.000162, -0.000077, 0, -0.000148, # cCA,cCB,cCC,cCD
60 -0.000167, -0.000162, -0.000148, 0, ], # cDA,cDB,cDC,cDD
61 )
62 crossTalkCoeffs2 = pexConfig.ListField(
63 dtype=float,
64 doc="crosstalk coefficients (secondary by the high-count pixel + 2pix)",
65 default=[0, 0.000051, 0.000050, 0.000053,
66 0.000051, 0, 0, 0.000050,
67 0.000050, 0, 0, 0.000051,
68 0.000053, 0.000050, 0.000051, 0, ],
69 )
71 relativeGainsPreampAndSigboard = pexConfig.ListField(
72 dtype=float,
73 doc="effective gain of combination of SIG board + preAmp for CCD, relatve to chB@CCD=5. \
74 g2*g3 in Yagi+2012",
75 default=[0.949, 0.993, 0.976, 0.996,
76 0.973, 0.984, 0.966, 0.977,
77 1.008, 0.989, 0.970, 0.976,
78 0.961, 0.966, 1.008, 0.967,
79 0.967, 0.984, 0.998, 1.000,
80 0.989, 1.000, 1.034, 1.030,
81 0.957, 1.019, 0.952, 0.979,
82 0.974, 1.015, 0.967, 0.962,
83 0.972, 0.932, 0.999, 0.963,
84 0.987, 0.985, 0.986, 1.012, ],
85 )
86 relativeGainsTotalBeforeOct2010 = pexConfig.ListField(
87 dtype=float,
88 doc="total effective gain of SIG board + preAmp + CCD on-chip amp, relatve to chB@CCD=5. \
89 effective before Oct 2010 g1*g2*g3 in Yagi+2012",
90 default=[1.0505, 1.06294, 1.08287, 1.06868,
91 1.06962, 1.13801, 1.21669, 1.21205,
92 0.99456, 1.01931, 1.038, 0.947461,
93 1.02245, 1.03952, 1.04192, 1.1544,
94 1.0055, 1.12767, 1.08451, 1.03309,
95 1.01532, 1.00000, 0.972106, 1.08237,
96 1.19014, 0.987111, 0.984164, 1.03382,
97 0.971049, 1.05928, 1.06713, 0.980967,
98 1.0107, 1.25209, 1.27565, 1.17183,
99 0.993771, 1.077, 1.03909, 1.03241, ],
100 )
102 relativeGainsTotalAfterOct2010 = pexConfig.ListField(
103 dtype=float,
104 doc="total effective gain of SIG board + preAmp + CCD on-chip amp, relatve to chB@CCD=5. \
105 effective after Oct 2010 g1*g2*g3 in Yagi+2012",
106 default=[1.0505, 1.06294, 1.08287, 1.06868,
107 1.06962, 1.13801, 1.21669, 1.21205,
108 0.99456, 1.01931, 1.038, 0.947461,
109 1.02245, 1.03952, 1.04192, 1.1544,
110 1.0055, 1.12767, 1.08451, 1.03309,
111 1.01532, 1.00000, 0.972106, 1.08237,
112 1.19014, 0.987111, 0.984164, 1.03382,
113 0.971049, 1.05928, 1.06713, 0.980967,
114 1.0107, 1.25209, 1.27565, 1.17183,
115 1.13003, 1.077, 1.03909, 1.03241, ],
116 )
118 shapeGainsArray = pexConfig.ListField(
119 dtype=int,
120 doc="Shape of gain arrays",
121 default=[10, 4], # 10 CCDs and 4 channels per CCD
122 minLength=1, # really 2, but there's a bug in pex_config
123 maxLength=2,
124 )
126 shapeCoeffsArray = pexConfig.ListField(
127 dtype=int,
128 doc="Shape of coeffs arrays",
129 default=[4, 4], # 4 channels and 4 mirror-symmetry patterns per channel
130 minLength=1, # really 2, but there's a bug in pex_config
131 maxLength=2,
132 )
134 def getCoeffs1(self):
135 """Return a 2-D numpy array of crosstalk coefficients of the proper
136 shape"""
137 return np.array(self.crossTalkCoeffs1).reshape(self.shapeCoeffsArray).tolist()
139 def getCoeffs2(self):
140 """Return a 2-D numpy array of crosstalk coefficients of the proper
141 shape"""
142 return np.array(self.crossTalkCoeffs2).reshape(self.shapeCoeffsArray).tolist()
144 def getGainsPreampSigboard(self):
145 """Return a 2-D numpy array of effective gain for preamp+SIG of the
146 proper shape"""
147 return np.array(self.relativeGainsPreampAndSigboard).reshape(self.shapeGainsArray).tolist()
149 def getGainsTotal(self, dateobs='2008-08-01'):
150 """Return a 2-D numpy array of effective total gain of the proper
151 shape"""
152 if dateobs < '2010-10': # may need rewritten to a more reliable way
153 return np.array(self.relativeGainsTotalBeforeOct2010).reshape(self.shapeGainsArray).tolist()
154 else:
155 return np.array(self.relativeGainsTotalAfterOct2010).reshape(self.shapeGainsArray).tolist()
158nAmp = 4
161def getCrosstalkX1(x, xmax=512):
162 """
163 Return the primary X positions in CCD which affect the count of input
164 x pixel
166 Parameters
167 ----------
168 x :
169 xpos in Ccd, affected by the returned pixels
170 xmax :
171 nx per amp
172 """
173 ctx1List = [x,
174 2 * xmax - x - 1,
175 2 * xmax + x,
176 4 * xmax - x - 1, ]
178 return ctx1List
181def getCrosstalkX2(x, xmax=512):
182 """
183 Return the 2ndary X positions (dX=2) in CCD which affect the count of
184 input x pixel
186 Those X pixels are read by 2-pixel ealier than the primary pixels.
188 Parameters
189 ----------
190 x :
191 xpos in Ccd, affected by the returned pixels
192 xmax :
193 nx per amp
194 """
195 # ch0,ch2: ctx2 = ctx1 - 2, ch1,ch3: ctx2 = ctx1 + 2
196 ctx2List = [x - 2,
197 2 * xmax - x - 1 + 2,
198 2 * xmax + x - 2,
199 4 * xmax - x - 1 + 2, ]
201 return ctx2List
204def getAmplifier(image, amp, ampReference=None, offset=2):
205 """Extract an image of the amplifier from the CCD, along with an offset
206 version
208 The amplifier image will be flipped (if required) to match the
209 orientation of a nominated reference amplifier.
211 An additional image, with the nominated offset applied, is also produced.
213 Parameters
214 ----------
215 image :
216 Image of CCD
217 amp :
218 Index of amplifier
219 ampReference :
220 Index of reference amplifier
221 offset :
222 Offset to apply
224 Returns
225 -------
226 amp_image :
227 amplifier image, offset amplifier image
228 """
229 height = image.getHeight()
230 ampBox = geom.Box2I(geom.Point2I(amp*512, 0), geom.Extent2I(512, height))
231 ampImage = image.Factory(image, ampBox, afwImage.LOCAL)
232 if ampReference is not None and amp % 2 != ampReference % 2:
233 ampImage = afwMath.flipImage(ampImage, True, False)
234 offBox = geom.Box2I(geom.Point2I(offset if amp == ampReference else 0, 0),
235 geom.Extent2I(510, height))
236 offsetImage = ampImage.Factory(ampImage, offBox, afwImage.LOCAL)
237 return ampImage, offsetImage
240def subtractCrosstalkYagi(mi, coeffs1List, coeffs2List, gainsPreampSig, minPixelToMask=45000,
241 crosstalkStr="CROSSTALK"):
242 """Subtract the crosstalk from MaskedImage mi given a set of coefficients
244 Based on procedure presented in Yagi et al. 2012, PASP in publication;
245 arXiv:1210.8212 The pixels affected by signal over minPixelToMask have
246 the crosstalkStr bit set
247 """
249 #
250 # These are the pixels that are bright enough to cause crosstalk (more
251 # precisely, the ones that we label as causing crosstalk; in reality all
252 # pixels cause crosstalk)
253 #
254 if True:
255 tempStr = "TEMP" # mask plane used to record the bright pixels that we need to mask
256 mi.getMask().addMaskPlane(tempStr)
257 fs = afwDetect.FootprintSet(mi, afwDetect.Threshold(minPixelToMask), tempStr)
259 mi.getMask().addMaskPlane(crosstalkStr)
260 afwDisplay.setMaskPlaneColor(crosstalkStr, afwDisplay.MAGENTA)
261 fs.setMask(mi.getMask(), crosstalkStr) # the crosstalkStr bit will now be set
262 # whenever we subtract crosstalk
263 crosstalk = mi.getMask().getPlaneBitMask(crosstalkStr)
265 if True:
266 # This python implementation is fairly fast
267 image = mi.getImage()
268 xtalk = afwImage.ImageF(image.getDimensions())
269 xtalk.set(0)
270 for i in range(4):
271 xAmp, xOff = getAmplifier(xtalk, i, i)
272 for j in range(4):
273 if i == j:
274 continue
275 gainRatio = gainsPreampSig[j] / gainsPreampSig[i]
276 jAmp, jOff = getAmplifier(image, j, i)
277 xAmp.scaledPlus(gainRatio*coeffs1List[i][j], jAmp)
278 xOff.scaledPlus(gainRatio*coeffs2List[i][j], jOff)
280 image -= xtalk
281 else:
282 nAmp = 4
283 subtractCrosstalk(mi, nAmp, coeffs1List, coeffs2List, gainsPreampSig)
285 #
286 # Clear the crosstalkStr bit in the original bright pixels, where tempStr
287 # is set
288 #
289 msk = mi.getMask()
290 temp = msk.getPlaneBitMask(tempStr)
291 xtalk_temp = crosstalk | temp
292 np_msk = msk.getArray()
293 np_msk[np.where(np.bitwise_and(np_msk, xtalk_temp) == xtalk_temp)] &= ~crosstalk
295 try:
296 msk.removeAndClearMaskPlane(tempStr, True) # added in afw #1853
297 except AttributeError:
298 afwDisplay.setMaskPlaneColor(tempStr, color="ignore")
301def main(visit=131634, ccd=None, threshold=45000, nSample=1, showCoeffs=True, fixXTalk=True,
302 plot=False, title=None):
303 if ccd is None:
304 visitList = list(range(nSample))
305 ccdList = ["simulated", ]
306 else:
307 ccdList = makeList(ccd)
308 visitList = makeList(visit)
310 coeffs = estimateCoeffs(visitList, ccdList, threshold=45000, plot=plot, title=title)
312 if showCoeffs:
313 printCoeffs(coeffs)
315 mi = readImage(visitList[0], ccdList[0])
316 if fixXTalk:
317 subtractXTalk(mi, coeffs, threshold)
319 return mi, coeffs
322class YagiCrosstalkConfig(pexConfig.Config):
323 coeffs = pexConfig.ConfigField(
324 dtype=CrosstalkYagiCoeffsConfig,
325 doc="Crosstalk coefficients by Yagi+ 2012",
326 )
327 crosstalkMaskPlane = pexConfig.Field(
328 dtype=str,
329 doc="Name of Mask plane for crosstalk corrected pixels",
330 default="CROSSTALK",
331 )
332 minPixelToMask = pexConfig.Field(
333 dtype=float,
334 doc="Minimum pixel value (in electrons) to cause crosstalkMaskPlane bit to be set",
335 default=45000,
336 )
339class YagiCrosstalkTask(pipeBase.Task):
340 ConfigClass = YagiCrosstalkConfig
342 def run(self, exposure):
343 coeffs1List = self.config.coeffs.getCoeffs1() # primary crosstalk
344 coeffs2List = self.config.coeffs.getCoeffs2() # secondary crosstalk
345 gainsPreampSig = self.config.coeffs.getGainsPreampSigboard()
346 if not np.any(coeffs1List):
347 self.log.info("No crosstalk info available. Skipping crosstalk corrections to CCD %s",
348 exposure.getDetector().getId())
349 return
351 self.log.info("Applying crosstalk corrections to CCD %s based on Yagi+2012",
352 exposure.getDetector().getId())
354 ccdId = int(exposure.getDetector().getId().getSerial())
355 gainsPreampSigCcd = gainsPreampSig[ccdId]
357 subtractCrosstalkYagi(exposure.getMaskedImage(), coeffs1List, coeffs2List,
358 gainsPreampSigCcd,
359 self.config.minPixelToMask, self.config.crosstalkMaskPlane)
362if __name__ == "__main__":
363 main()