LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
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, 2009, 2010 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 from __future__ import division
25 
26 from optparse import OptionParser
27 import lsst.afw.image as afwImage
28 import lsst.afw.geom as afwGeom
29 import lsst.afw.detection as afwDet
30 import lsst.pex.policy as pexPolicy
31 import lsst.pex.logging as pexLog
32 import lsst.daf.persistence as dafPersist
33 import lsst.daf.base as dafBase
34 import numpy as num
35 import lsst.afw.display.ds9 as ds9
36 import lsst.pex.config as pexConfig
37 
38 scaling = 5
39 
41  """Parse the command line options."""
42  parser = OptionParser(
43  usage="""%prog cdDiffSources crDiffExposure
44 
45 Read in sources and test for junk""")
46  options, args = parser.parse_args()
47  if len(args) != 2:
48  parser.error("incorrect number of arguments")
49  return options, args
50 
51 def readSourceSet(boostFile):
52  loc = dafPersist.LogicalLocation(boostFile)
53  storageList = dafPersist.StorageList()
54  additionalData = dafBase.PropertySet()
55  persistence = dafPersist.Persistence.getPersistence(pexPolicy.Policy())
56  storageList.append(persistence.getRetrieveStorage("BoostStorage", loc))
57  psvptr = persistence.unsafeRetrieve("PersistableSourceVector", storageList, additionalData)
58  psv = afwDet.PersistableSourceVector.swigConvert(psvptr)
59  return psv.getSources()
60 
61 class DiaSourceAnalystConfig(pexConfig.Config):
62  srcBadMaskPlanes = pexConfig.ListField(
63  dtype = str,
64  doc = """Mask planes that lead to an invalid detection.
65  Options: NO_DATA EDGE SAT BAD CR INTRP
66  E.g. : NO_DATA SAT BAD allows CR-masked and interpolated pixels""",
67  default = ("NO_DATA", "EDGE", "SAT", "BAD")
68  )
69  fBadPixels = pexConfig.Field(
70  dtype = float,
71  doc = "Fraction of bad pixels allowed in footprint",
72  default = 0.1
73  )
74  fluxPolarityRatio = pexConfig.Field(
75  dtype = float,
76  doc = "Minimum fraction of flux in correct-polarity pixels",
77  default = 0.75
78  )
79  nPolarityRatio = pexConfig.Field(
80  dtype = float,
81  doc = "Minimum fraction of correct-polarity pixels in unmasked subset",
82  default = 0.7
83  )
84  nMaskedRatio = pexConfig.Field(
85  dtype = float,
86  doc = "Minimum fraction of correct-polarity unmasked to masked pixels",
87  default = 0.6,
88  )
89  nGoodRatio = pexConfig.Field(
90  dtype = float,
91  doc = "Minimum fraction of correct-polarity unmasked to all pixels",
92  default = 0.5
93  )
94 
95 
96 class DiaSourceAnalyst(object):
97  def __init__(self, config):
98  self.config = config
99 
100  self.bitMask = 0
101  srcBadMaskPlanes = self.config.srcBadMaskPlanes
102  for maskPlane in srcBadMaskPlanes:
103  self.bitMask |= afwImage.MaskU_getPlaneBitMask(maskPlane)
104 
105  def countDetected(self, mask):
106  idxP = num.where(mask & afwImage.MaskU_getPlaneBitMask("DETECTED"))
107  idxN = num.where(mask & afwImage.MaskU_getPlaneBitMask("DETECTED_NEGATIVE"))
108  return len(idxP[0]), len(idxN[0])
109 
110  def countMasked(self, mask):
111  idxM = num.where(mask & self.bitMask)
112  return len(idxM[0])
113 
114  def countPolarity(self, mask, pixels):
115  unmasked = ((mask & self.bitMask) == 0)
116  idxP = num.where( (pixels >= 0) & unmasked)
117  idxN = num.where( (pixels < 0) & unmasked)
118  fluxP = num.sum(pixels[idxP])
119  fluxN = num.sum(pixels[idxN])
120  #import pdb; pdb.set_trace()
121  return len(idxP[0]), len(idxN[0]), fluxP, fluxN
122 
123  def testSource(self, source, subMi):
124  imArr, maArr, varArr = subMi.getArrays()
125  flux = source.getApFlux()
126 
127  nPixels = subMi.getWidth() * subMi.getHeight()
128  nPos, nNeg, fPos, fNeg = self.countPolarity(maArr, imArr)
129  nDetPos, nDetNeg = self.countDetected(maArr)
130  nMasked = self.countMasked(maArr)
131  assert (nPixels == (nMasked + nPos + nNeg))
132 
133  # 1) Too many pixels in the detection are masked
134  fMasked = (nMasked / nPixels)
135  fMaskedTol = self.config.fBadPixels
136  if fMasked > fMaskedTol:
137  pexLog.Trace("lsst.ip.diffim.DiaSourceAnalysis", 1,
138  "Candidate %d : BAD fBadPixels %.2f > %.2f" % (source.getId(), fMasked, fMaskedTol))
139  return False
140 
141  if flux > 0:
142  # positive-going source
143  fluxRatio = fPos / (fPos + abs(fNeg))
144  ngoodRatio = nPos / nPixels
145  maskRatio = nPos / (nPos + nMasked)
146  npolRatio = nPos / (nPos + nNeg)
147  else:
148  # negative-going source
149  fluxRatio = abs(fNeg)/ (fPos + abs(fNeg))
150  ngoodRatio = nNeg / nPixels
151  maskRatio = nNeg / (nNeg + nMasked)
152  npolRatio = nNeg / (nNeg + nPos)
153 
154  # 2) Not enough flux in unmasked correct-polarity pixels
155  fluxRatioTolerance = self.config.fluxPolarityRatio
156  if fluxRatio < fluxRatioTolerance:
157  pexLog.Trace("lsst.ip.diffim.DiaSourceAnalysis", 1,
158  "Candidate %d : BAD flux polarity %.2f < %.2f (pos=%.2f neg=%.2f)" % (source.getId(),
159  fluxRatio, fluxRatioTolerance, fPos, fNeg))
160  return False
161 
162  # 3) Not enough unmasked pixels of correct polarity
163  polarityTolerance = self.config.nPolarityRatio
164  if npolRatio < polarityTolerance:
165  pexLog.Trace("lsst.ip.diffim.DiaSourceAnalysis", 1,
166  "Candidate %d : BAD polarity count %.2f < %.2f (pos=%d neg=%d)" % (source.getId(),
167  npolRatio, polarityTolerance, nPos, nNeg))
168  return False
169 
170  # 4) Too many masked vs. correct polarity pixels
171  maskedTolerance = self.config.nMaskedRatio
172  if maskRatio < maskedTolerance:
173  pexLog.Trace("lsst.ip.diffim.DiaSourceAnalysis", 1,
174  "Candidate %d : BAD unmasked count %.2f < %.2f (pos=%d neg=%d mask=%d)" % (source.getId(),
175  maskRatio, maskedTolerance, nPos, nNeg, nMasked))
176  return False
177 
178  # 5) Too few unmasked, correct polarity pixels
179  ngoodTolerance = self.config.nGoodRatio
180  if ngoodRatio < ngoodTolerance:
181  pexLog.Trace("lsst.ip.diffim.DiaSourceAnalysis", 1,
182  "Candidate %d : BAD good pixel count %.2f < %.2f (pos=%d neg=%d tot=%d)" % (source.getId(),
183  ngoodRatio, ngoodTolerance, nPos, nNeg, nPixels))
184  return False
185 
186  pexLog.Trace("lsst.ip.diffim.DiaSourceAnalysis", 1,
187  "Candidate %d : OK flux=%.2f nPos=%d nNeg=%d nTot=%d nDetPos=%d nDetNeg=%d fPos=%.2f fNeg=%2f" % (source.getId(),
188  flux, nPos, nNeg, nPixels, nDetPos, nDetNeg, fPos, fNeg))
189  return True
190 
191 
192 def main():
193  """Main program"""
194  options, args = parseOptions()
195  (crDiffSourceFile, crDiffExposureFile) = args
196 
197  crDiffSources = readSourceSet(crDiffSourceFile)
198  crDiffExposure = afwImage.ExposureF(crDiffExposureFile)
199  #import pdb; pdb.set_trace()
200 
201  analyst = DiaSourceAnalyst()
202 
203  expX0 = crDiffExposure.getX0()
204  expY0 = crDiffExposure.getY0()
205  expX1 = expX0 + crDiffExposure.getWidth() - 1
206  expY1 = expY0 + crDiffExposure.getHeight() - 1
207 
208  for i in range(crDiffSources.size()):
209  crDiffSource = crDiffSources[i]
210 
211  # This segfaults; revisit once the stack stabilizes
212  #footprint = crDiffSource.getFootprint()
213  #bbox = footprint.getBBox()
214 
215  xAstrom = crDiffSource.getXAstrom()
216  yAstrom = crDiffSource.getYAstrom()
217  Ixx = max(1.0, crDiffSource.getIxx())
218  Iyy = max(1.0, crDiffSource.getIyy())
219  x0 = max(expX0, int(xAstrom - scaling * Ixx))
220  x1 = min(expX1, int(xAstrom + scaling * Ixx))
221  y0 = max(expY0, int(yAstrom - scaling * Iyy))
222  y1 = min(expY1, int(yAstrom + scaling * Iyy))
223  bbox = afwGeom.Box2I(afwGeom.Point2I(x0, y0),
224  afwGeom.Point2I(x1, y1))
225  subExp = afwImage.ExposureF(crDiffExposure, bbox)
226  subMi = subExp.getMaskedImage()
227  imArr, maArr, varArr = subMi.getArrays()
228 
229  if analyst.testSource(crDiffSource, subMi):
230  ds9.mtv(subExp, frame=1)
231  raw_input('Next: ')
232 if __name__ == "__main__":
233  main()
Class for logical location of a persisted Persistable instance.
a container for holding hierarchical configuration data in memory.
Definition: Policy.h:169
limited backward compatibility to the DC2 run-time trace facilities
Definition: Trace.h:93
An integer coordinate rectangle.
Definition: Box.h:53
Class for storing generic metadata.
Definition: PropertySet.h:82