LSSTApplications  18.1.0
LSSTDataManagementBasePackage
diaSourceAnalysis.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 #
4 # LSST Data Management System
5 # Copyright 2008-2016 LSST Corporation.
6 #
7 # This product includes software developed by the
8 # LSST Project (http://www.lsst.org/).
9 #
10 # This program is free software: you can redistribute it and/or modify
11 # it under the terms of the GNU General Public License as published by
12 # the Free Software Foundation, either version 3 of the License, or
13 # (at your option) any later version.
14 #
15 # This program is distributed in the hope that it will be useful,
16 # but WITHOUT ANY WARRANTY; without even the implied warranty of
17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 # GNU General Public License for more details.
19 #
20 # You should have received a copy of the LSST License Statement and
21 # the GNU General Public License along with this program. If not,
22 # see <http://www.lsstcorp.org/LegalNotices/>.
23 #
24 
25 __all__ = ["DiaSourceAnalystConfig", "DiaSourceAnalyst"]
26 
27 import lsst.afw.image as afwImage
28 from lsst.log import Log
29 import numpy as num
30 import lsst.pex.config as pexConfig
31 
32 scaling = 5
33 
34 
35 class DiaSourceAnalystConfig(pexConfig.Config):
36  srcBadMaskPlanes = pexConfig.ListField(
37  dtype=str,
38  doc="""Mask planes that lead to an invalid detection.
39  Options: NO_DATA EDGE SAT BAD CR INTRP
40  E.g. : NO_DATA SAT BAD allows CR-masked and interpolated pixels""",
41  default=("NO_DATA", "EDGE", "SAT", "BAD")
42  )
43  fBadPixels = pexConfig.Field(
44  dtype=float,
45  doc="Fraction of bad pixels allowed in footprint",
46  default=0.1
47  )
48  fluxPolarityRatio = pexConfig.Field(
49  dtype=float,
50  doc="Minimum fraction of flux in correct-polarity pixels",
51  default=0.75
52  )
53  nPolarityRatio = pexConfig.Field(
54  dtype=float,
55  doc="Minimum fraction of correct-polarity pixels in unmasked subset",
56  default=0.7
57  )
58  nMaskedRatio = pexConfig.Field(
59  dtype=float,
60  doc="Minimum fraction of correct-polarity unmasked to masked pixels",
61  default=0.6,
62  )
63  nGoodRatio = pexConfig.Field(
64  dtype=float,
65  doc="Minimum fraction of correct-polarity unmasked to all pixels",
66  default=0.5
67  )
68 
69 
71 
72  def __init__(self, config):
73  self.config = config
74  self.log = Log.getLogger("ip.diffim.DiaSourceAnalysis")
75 
76  self.bitMask = 0
77  srcBadMaskPlanes = self.config.srcBadMaskPlanes
78  for maskPlane in srcBadMaskPlanes:
79  self.bitMask |= afwImage.Mask.getPlaneBitMask(maskPlane)
80 
81  def countDetected(self, mask):
82  idxP = num.where(mask & afwImage.Mask.getPlaneBitMask("DETECTED"))
83  idxN = num.where(mask & afwImage.Mask.getPlaneBitMask("DETECTED_NEGATIVE"))
84  return len(idxP[0]), len(idxN[0])
85 
86  def countMasked(self, mask):
87  idxM = num.where(mask & self.bitMask)
88  return len(idxM[0])
89 
90  def countPolarity(self, mask, pixels):
91  unmasked = ((mask & self.bitMask) == 0)
92  idxP = num.where((pixels >= 0) & unmasked)
93  idxN = num.where((pixels < 0) & unmasked)
94  fluxP = num.sum(pixels[idxP])
95  fluxN = num.sum(pixels[idxN])
96  return len(idxP[0]), len(idxN[0]), fluxP, fluxN
97 
98  def testSource(self, source, subMi):
99  imArr, maArr, varArr = subMi.getArrays()
100  flux = source.getApFlux()
101 
102  nPixels = subMi.getWidth() * subMi.getHeight()
103  nPos, nNeg, fPos, fNeg = self.countPolarity(maArr, imArr)
104  nDetPos, nDetNeg = self.countDetected(maArr)
105  nMasked = self.countMasked(maArr)
106  assert(nPixels == (nMasked + nPos + nNeg))
107 
108  # 1) Too many pixels in the detection are masked
109  fMasked = (nMasked / nPixels)
110  fMaskedTol = self.config.fBadPixels
111  if fMasked > fMaskedTol:
112  self.log.debug("Candidate %d : BAD fBadPixels %.2f > %.2f", source.getId(), fMasked, fMaskedTol)
113  return False
114 
115  if flux > 0:
116  # positive-going source
117  fluxRatio = fPos / (fPos + abs(fNeg))
118  ngoodRatio = nPos / nPixels
119  maskRatio = nPos / (nPos + nMasked)
120  npolRatio = nPos / (nPos + nNeg)
121  else:
122  # negative-going source
123  fluxRatio = abs(fNeg) / (fPos + abs(fNeg))
124  ngoodRatio = nNeg / nPixels
125  maskRatio = nNeg / (nNeg + nMasked)
126  npolRatio = nNeg / (nNeg + nPos)
127 
128  # 2) Not enough flux in unmasked correct-polarity pixels
129  fluxRatioTolerance = self.config.fluxPolarityRatio
130  if fluxRatio < fluxRatioTolerance:
131  self.log.debug("Candidate %d : BAD flux polarity %.2f < %.2f (pos=%.2f neg=%.2f)",
132  source.getId(), fluxRatio, fluxRatioTolerance, fPos, fNeg)
133  return False
134 
135  # 3) Not enough unmasked pixels of correct polarity
136  polarityTolerance = self.config.nPolarityRatio
137  if npolRatio < polarityTolerance:
138  self.log.debug("Candidate %d : BAD polarity count %.2f < %.2f (pos=%d neg=%d)",
139  source.getId(), npolRatio, polarityTolerance, nPos, nNeg)
140  return False
141 
142  # 4) Too many masked vs. correct polarity pixels
143  maskedTolerance = self.config.nMaskedRatio
144  if maskRatio < maskedTolerance:
145  self.log.debug("Candidate %d : BAD unmasked count %.2f < %.2f (pos=%d neg=%d mask=%d)",
146  source.getId(), maskRatio, maskedTolerance, nPos, nNeg, nMasked)
147  return False
148 
149  # 5) Too few unmasked, correct polarity pixels
150  ngoodTolerance = self.config.nGoodRatio
151  if ngoodRatio < ngoodTolerance:
152  self.log.debug("Candidate %d : BAD good pixel count %.2f < %.2f (pos=%d neg=%d tot=%d)",
153  source.getId(), ngoodRatio, ngoodTolerance, nPos, nNeg, nPixels)
154  return False
155 
156  self.log.debug("Candidate %d : OK flux=%.2f nPos=%d nNeg=%d nTot=%d nDetPos=%d nDetNeg=%d "
157  "fPos=%.2f fNeg=%2f",
158  source.getId(), flux, nPos, nNeg, nPixels, nDetPos, nDetNeg, fPos, fNeg)
159  return True
Angle abs(Angle const &a)
Definition: Angle.h:106
Definition: Log.h:691
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...