LSSTApplications  18.1.0
LSSTDataManagementBasePackage
defects.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 #
4 # Copyright 2008-2017 AURA/LSST.
5 #
6 # This product includes software developed by the
7 # LSST Project (http://www.lsst.org/).
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 LSST License Statement and
20 # the GNU General Public License along with this program. If not,
21 # see <https://www.lsstcorp.org/LegalNotices/>.
22 #
23 """Support for image defects"""
24 
25 __all__ = ("Defects",)
26 
27 import logging
28 import itertools
29 import collections.abc
30 from deprecated.sphinx import deprecated
31 import numpy as np
32 import copy
33 import datetime
34 import math
35 import numbers
36 import os.path
37 import astropy.table
38 
39 import lsst.geom
40 import lsst.pex.policy as policy
41 import lsst.afw.table
42 import lsst.afw.detection
43 import lsst.afw.image
44 import lsst.afw.geom
45 from lsst.daf.base import PropertyList
46 
47 from . import Defect
48 
49 log = logging.getLogger(__name__)
50 
51 SCHEMA_NAME_KEY = "DEFECTS_SCHEMA"
52 SCHEMA_VERSION_KEY = "DEFECTS_SCHEMA_VERSION"
53 
54 
55 @deprecated(reason="Policy defect files no longer supported (will be removed after v18)",
56  category=FutureWarning)
57 def policyToBadRegionList(policyFile):
58  """Given a Policy file describing a CCD's bad pixels, return a vector of BadRegion::Ptr"""
59 
60  badPixelsPolicy = policy.Policy.createPolicy(policyFile)
61  badPixels = []
62 
63  if badPixelsPolicy.exists("Defects"):
64  d = badPixelsPolicy.getArray("Defects")
65  for reg in d:
66  x0 = reg.get("x0")
67  width = reg.get("width")
68  if not width:
69  x1 = reg.get("x1")
70  width = x1 - x0 - 1
71 
72  y0 = reg.get("y0")
73  if reg.exists("height"):
74  height = reg.get("height")
75  else:
76  y1 = reg.get("y1")
77  height = y1 - y0 - 1
78 
79  bbox = lsst.geom.BoxI(lsst.geom.PointI(x0, y0), lsst.geom.ExtentI(width, height))
80  badPixels.append(Defect(bbox))
81 
82  return badPixels
83 
84 
85 class Defects(collections.abc.MutableSequence):
86  """Collection of `lsst.meas.algorithms.Defect`.
87 
88  Parameters
89  ----------
90  defectList : iterable of `lsst.meas.algorithms.Defect`
91  or `lsst.geom.BoxI`, optional
92  Collections of defects to apply to the image.
93  """
94 
95  _OBSTYPE = "defects"
96  """The calibration type used for ingest."""
97 
98  def __init__(self, defectList=None, metadata=None):
99  self._defects = []
100 
101  if metadata is not None:
102  self._metadata = metadata
103  else:
104  self.setMetadata()
105 
106  if defectList is None:
107  return
108 
109  # Ensure that type checking
110  for d in defectList:
111  self.append(d)
112 
113  def _check_value(self, value):
114  """Check that the supplied value is a `~lsst.meas.algorithms.Defect`
115  or can be converted to one.
116 
117  Parameters
118  ----------
119  value : `object`
120  Value to check.
121 
122  Returns
123  -------
124  new : `~lsst.meas.algorithms.Defect`
125  Either the supplied value or a new object derived from it.
126 
127  Raises
128  ------
129  ValueError
130  Raised if the supplied value can not be converted to
131  `~lsst.meas.algorithms.Defect`
132  """
133  if isinstance(value, Defect):
134  pass
135  elif isinstance(value, lsst.geom.BoxI):
136  value = Defect(value)
137  elif isinstance(value, lsst.geom.PointI):
138  value = Defect(lsst.geom.Box2I(value, lsst.geom.Extent2I(1, 1)))
139  elif isinstance(value, lsst.afw.image.DefectBase):
140  value = Defect(value.getBBox())
141  else:
142  raise ValueError(f"Defects must be of type Defect, BoxI, or PointI, not '{value!r}'")
143  return value
144 
145  def __len__(self):
146  return len(self._defects)
147 
148  def __getitem__(self, index):
149  return self._defects[index]
150 
151  def __setitem__(self, index, value):
152  """Can be given a `~lsst.meas.algorithms.Defect` or a `lsst.geom.BoxI`
153  """
154  self._defects[index] = self._check_value(value)
155 
156  def __iter__(self):
157  return iter(self._defects)
158 
159  def __delitem__(self, index):
160  del self._defects[index]
161 
162  def __eq__(self, other):
163  """Compare if two `Defects` are equal.
164 
165  Two `Defects` are equal if their bounding boxes are equal and in
166  the same order. Metadata content is ignored.
167  """
168  if not isinstance(other, self.__class__):
169  return False
170 
171  # checking the bboxes with zip() only works if same length
172  if len(self) != len(other):
173  return False
174 
175  # Assume equal if bounding boxes are equal
176  for d1, d2 in zip(self, other):
177  if d1.getBBox() != d2.getBBox():
178  return False
179 
180  return True
181 
182  def __str__(self):
183  return "Defects(" + ",".join(str(d.getBBox()) for d in self) + ")"
184 
185  def insert(self, index, value):
186  self._defects.insert(index, self._check_value(value))
187 
188  def getMetadata(self):
189  """Retrieve metadata associated with these `Defects`.
190 
191  Returns
192  -------
193  meta : `lsst.daf.base.PropertyList`
194  Metadata. The returned `~lsst.daf.base.PropertyList` can be
195  modified by the caller and the changes will be written to
196  external files.
197  """
198  return self._metadata
199 
200  def setMetadata(self, metadata=None):
201  """Store a copy of the supplied metadata with the defects.
202 
203  Parameters
204  ----------
205  metadata : `lsst.daf.base.PropertyList`, optional
206  Metadata to associate with the defects. Will be copied and
207  overwrite existing metadata. If not supplied the existing
208  metadata will be reset.
209  """
210  if metadata is None:
211  self._metadata = PropertyList()
212  else:
213  self._metadata = copy.copy(metadata)
214 
215  # Ensure that we have the obs type required by calibration ingest
216  self._metadata["OBSTYPE"] = self._OBSTYPE
217 
218  def copy(self):
219  """Copy the defects to a new list, creating new defects from the
220  bounding boxes.
221 
222  Returns
223  -------
224  new : `Defects`
225  New list with new `Defect` entries.
226 
227  Notes
228  -----
229  This is not a shallow copy in that new `Defect` instances are
230  created from the original bounding boxes. It's also not a deep
231  copy since the bounding boxes are not recreated.
232  """
233  return self.__class__(d.getBBox() for d in self)
234 
235  def transpose(self):
236  """Make a transposed copy of this defect list.
237 
238  Returns
239  -------
240  retDefectList : `Defects`
241  Transposed list of defects.
242  """
243  retDefectList = self.__class__()
244  for defect in self:
245  bbox = defect.getBBox()
246  dimensions = bbox.getDimensions()
247  nbbox = lsst.geom.Box2I(lsst.geom.Point2I(bbox.getMinY(), bbox.getMinX()),
248  lsst.geom.Extent2I(dimensions[1], dimensions[0]))
249  retDefectList.append(nbbox)
250  return retDefectList
251 
252  def maskPixels(self, maskedImage, maskName="BAD"):
253  """Set mask plane based on these defects.
254 
255  Parameters
256  ----------
257  maskedImage : `lsst.afw.image.MaskedImage`
258  Image to process. Only the mask plane is updated.
259  maskName : str, optional
260  Mask plane name to use.
261  """
262  # mask bad pixels
263  mask = maskedImage.getMask()
264  bitmask = mask.getPlaneBitMask(maskName)
265  for defect in self:
266  bbox = defect.getBBox()
267  lsst.afw.geom.SpanSet(bbox).clippedTo(mask.getBBox()).setMask(mask, bitmask)
268 
269  def toFitsRegionTable(self):
270  """Convert defect list to `~lsst.afw.table.BaseCatalog` using the
271  FITS region standard.
272 
273  Returns
274  -------
275  table : `lsst.afw.table.BaseCatalog`
276  Defects in tabular form.
277 
278  Notes
279  -----
280  The table created uses the
281  `FITS regions <https://fits.gsfc.nasa.gov/registry/region.html>`_
282  definition tabular format. The ``X`` and ``Y`` coordinates are
283  converted to FITS Physical coordinates that have origin pixel (1, 1)
284  rather than the (0, 0) used in LSST software.
285  """
286  schema = lsst.afw.table.Schema()
287  x = schema.addField("X", type="D", units="pix", doc="X coordinate of center of shape")
288  y = schema.addField("Y", type="D", units="pix", doc="Y coordinate of center of shape")
289  shape = schema.addField("SHAPE", type="String", size=16, doc="Shape defined by these values")
290  r = schema.addField("R", type="ArrayD", size=2, units="pix", doc="Extents")
291  rotang = schema.addField("ROTANG", type="D", units="deg", doc="Rotation angle")
292  component = schema.addField("COMPONENT", type="I", doc="Index of this region")
293  table = lsst.afw.table.BaseCatalog(schema)
294  table.resize(len(self._defects))
295 
296  for i, defect in enumerate(self._defects):
297  box = defect.getBBox()
298  # Correct for the FITS 1-based offset
299  table[i][x] = box.getCenterX() + 1.0
300  table[i][y] = box.getCenterY() + 1.0
301  width = box.getWidth()
302  height = box.getHeight()
303 
304  if width == 1 and height == 1:
305  # Call this a point
306  shapeType = "POINT"
307  else:
308  shapeType = "BOX"
309  table[i][shape] = shapeType
310  table[i][r] = np.array([width, height], dtype=np.float64)
311  table[i][rotang] = 0.0
312  table[i][component] = i
313 
314  # Set some metadata in the table (force OBSTYPE to exist)
315  metadata = copy.copy(self.getMetadata())
316  metadata["OBSTYPE"] = self._OBSTYPE
317  metadata[SCHEMA_NAME_KEY] = "FITS Region"
318  metadata[SCHEMA_VERSION_KEY] = 1
319  table.setMetadata(metadata)
320 
321  return table
322 
323  def writeFits(self, *args):
324  """Write defect list to FITS.
325 
326  Parameters
327  ----------
328  *args
329  Arguments to be forwarded to
330  `lsst.afw.table.BaseCatalog.writeFits`.
331  """
332  table = self.toFitsRegionTable()
333 
334  # Add some additional headers useful for tracking purposes
335  metadata = table.getMetadata()
336  now = datetime.datetime.utcnow()
337  metadata["DATE"] = now.isoformat()
338  metadata["CALIB_CREATION_DATE"] = now.strftime("%Y-%m-%d")
339  metadata["CALIB_CREATION_TIME"] = now.strftime("%T %Z").strip()
340 
341  table.writeFits(*args)
342 
343  def toSimpleTable(self):
344  """Convert defects to a simple table form that we use to write
345  to text files.
346 
347  Returns
348  -------
349  table : `lsst.afw.table.BaseCatalog`
350  Defects in simple tabular form.
351 
352  Notes
353  -----
354  These defect tables are used as the human readable definitions
355  of defects in calibration data definition repositories. The format
356  is to use four columns defined as follows:
357 
358  x0 : `int`
359  X coordinate of bottom left corner of box.
360  y0 : `int`
361  Y coordinate of bottom left corner of box.
362  width : `int`
363  X extent of the box.
364  height : `int`
365  Y extent of the box.
366  """
367  schema = lsst.afw.table.Schema()
368  x = schema.addField("x0", type="I", units="pix",
369  doc="X coordinate of bottom left corner of box")
370  y = schema.addField("y0", type="I", units="pix",
371  doc="Y coordinate of bottom left corner of box")
372  width = schema.addField("width", type="I", units="pix",
373  doc="X extent of box")
374  height = schema.addField("height", type="I", units="pix",
375  doc="Y extent of box")
376  table = lsst.afw.table.BaseCatalog(schema)
377  table.resize(len(self._defects))
378 
379  for i, defect in enumerate(self._defects):
380  box = defect.getBBox()
381  table[i][x] = box.getBeginX()
382  table[i][y] = box.getBeginY()
383  table[i][width] = box.getWidth()
384  table[i][height] = box.getHeight()
385 
386  # Set some metadata in the table (force OBSTYPE to exist)
387  metadata = copy.copy(self.getMetadata())
388  metadata["OBSTYPE"] = self._OBSTYPE
389  metadata[SCHEMA_NAME_KEY] = "Simple"
390  metadata[SCHEMA_VERSION_KEY] = 1
391  table.setMetadata(metadata)
392 
393  return table
394 
395  def writeText(self, filename):
396  """Write the defects out to a text file with the specified name.
397 
398  Parameters
399  ----------
400  filename : `str`
401  Name of the file to write. The file extension ".ecsv" will
402  always be used.
403 
404  Returns
405  -------
406  used : `str`
407  The name of the file used to write the data (which may be
408  different from the supplied name given the change to file
409  extension).
410 
411  Notes
412  -----
413  The file is written to ECSV format and will include any metadata
414  associated with the `Defects`.
415  """
416 
417  # Using astropy table is the easiest way to serialize to ecsv
418  afwTable = self.toSimpleTable()
419  table = afwTable.asAstropy()
420 
421  metadata = afwTable.getMetadata()
422  now = datetime.datetime.utcnow()
423  metadata["DATE"] = now.isoformat()
424  metadata["CALIB_CREATION_DATE"] = now.strftime("%Y-%m-%d")
425  metadata["CALIB_CREATION_TIME"] = now.strftime("%T %Z").strip()
426 
427  table.meta = metadata.toDict()
428 
429  # Force file extension to .ecsv
430  path, ext = os.path.splitext(filename)
431  filename = path + ".ecsv"
432  table.write(filename, format="ascii.ecsv")
433  return filename
434 
435  @staticmethod
436  def _get_values(values, n=1):
437  """Retrieve N values from the supplied values.
438 
439  Parameters
440  ----------
441  values : `numbers.Number` or `list` or `np.array`
442  Input values.
443  n : `int`
444  Number of values to retrieve.
445 
446  Returns
447  -------
448  vals : `list` or `np.array` or `numbers.Number`
449  Single value from supplied list if ``n`` is 1, or `list`
450  containing first ``n`` values from supplied values.
451 
452  Notes
453  -----
454  Some supplied tables have vectors in some columns that can also
455  be scalars. This method can be used to get the first number as
456  a scalar or the first N items from a vector as a vector.
457  """
458  if n == 1:
459  if isinstance(values, numbers.Number):
460  return values
461  else:
462  return values[0]
463 
464  return values[:n]
465 
466  @classmethod
467  def fromTable(cls, table):
468  """Construct a `Defects` from the contents of a
469  `~lsst.afw.table.BaseCatalog`.
470 
471  Parameters
472  ----------
473  table : `lsst.afw.table.BaseCatalog`
474  Table with one row per defect.
475 
476  Returns
477  -------
478  defects : `Defects`
479  A `Defects` list.
480 
481  Notes
482  -----
483  Two table formats are recognized. The first is the
484  `FITS regions <https://fits.gsfc.nasa.gov/registry/region.html>`_
485  definition tabular format written by `toFitsRegionTable` where the
486  pixel origin is corrected from FITS 1-based to a 0-based origin.
487  The second is the legacy defects format using columns ``x0``, ``y0``
488  (bottom left hand pixel of box in 0-based coordinates), ``width``
489  and ``height``.
490 
491  The FITS standard regions can only read BOX, POINT, or ROTBOX with
492  a zero degree rotation.
493  """
494 
495  defectList = []
496 
497  schema = table.getSchema()
498 
499  # Check schema to see which definitions we have
500  if "X" in schema and "Y" in schema and "R" in schema and "SHAPE" in schema:
501  # This is a FITS region style table
502  isFitsRegion = True
503 
504  elif "x0" in schema and "y0" in schema and "width" in schema and "height" in schema:
505  # This is a classic LSST-style defect table
506  isFitsRegion = False
507 
508  else:
509  raise ValueError("Unsupported schema for defects extraction")
510 
511  for r in table:
512  record = r.extract("*")
513 
514  if isFitsRegion:
515  # Coordinates can be arrays (some shapes in the standard
516  # require this)
517  # Correct for FITS 1-based origin
518  xcen = cls._get_values(record["X"]) - 1.0
519  ycen = cls._get_values(record["Y"]) - 1.0
520  shape = record["SHAPE"].upper()
521  if shape == "BOX":
523  lsst.geom.Extent2I(cls._get_values(record["R"],
524  n=2)))
525  elif shape == "POINT":
526  # Handle the case where we have an externally created
527  # FITS file.
528  box = lsst.geom.Point2I(xcen, ycen)
529  elif shape == "ROTBOX":
530  # Astropy regions always writes ROTBOX
531  rotang = cls._get_values(record["ROTANG"])
532  # We can support 0 or 90 deg
533  if math.isclose(rotang % 90.0, 0.0):
534  # Two values required
535  r = cls._get_values(record["R"], n=2)
536  if math.isclose(rotang % 180.0, 0.0):
537  width = r[0]
538  height = r[1]
539  else:
540  width = r[1]
541  height = r[0]
543  lsst.geom.Extent2I(width, height))
544  else:
545  log.warning("Defect can not be defined using ROTBOX with non-aligned rotation angle")
546  continue
547  else:
548  log.warning("Defect lists can only be defined using BOX or POINT not %s", shape)
549  continue
550 
551  elif "x0" in record and "y0" in record and "width" in record and "height" in record:
552  # This is a classic LSST-style defect table
553  box = lsst.geom.Box2I(lsst.geom.Point2I(record["x0"], record["y0"]),
554  lsst.geom.Extent2I(record["width"], record["height"]))
555 
556  defectList.append(box)
557 
558  defects = cls(defectList)
559  defects.setMetadata(table.getMetadata())
560 
561  # Once read, the schema headers are irrelevant
562  metadata = defects.getMetadata()
563  for k in (SCHEMA_NAME_KEY, SCHEMA_VERSION_KEY):
564  if k in metadata:
565  del metadata[k]
566 
567  return defects
568 
569  @classmethod
570  def readFits(cls, *args):
571  """Read defect list from FITS table.
572 
573  Parameters
574  ----------
575  *args
576  Arguments to be forwarded to
577  `lsst.afw.table.BaseCatalog.writeFits`.
578 
579  Returns
580  -------
581  defects : `Defects`
582  Defects read from a FITS table.
583  """
584  table = lsst.afw.table.BaseCatalog.readFits(*args)
585  return cls.fromTable(table)
586 
587  @classmethod
588  def readText(cls, filename):
589  """Read defect list from standard format text table file.
590 
591  Parameters
592  ----------
593  filename : `str`
594  Name of the file containing the defects definitions.
595 
596  Returns
597  -------
598  defects : `Defects`
599  Defects read from a FITS table.
600  """
601  table = astropy.table.Table.read(filename)
602 
603  # Need to convert the Astropy table to afw table
604  schema = lsst.afw.table.Schema()
605  for colName in table.columns:
606  schema.addField(colName, units=str(table[colName].unit),
607  type=table[colName].dtype.type)
608 
609  # Create AFW table that is required by fromTable()
610  afwTable = lsst.afw.table.BaseCatalog(schema)
611 
612  afwTable.resize(len(table))
613  for colName in table.columns:
614  # String columns will fail -- currently we do not expect any
615  afwTable[colName] = table[colName]
616 
617  # Copy in the metadata from the astropy table
618  metadata = PropertyList()
619  for k, v in table.meta.items():
620  metadata[k] = v
621  afwTable.setMetadata(metadata)
622 
623  # Extract defect information from the table itself
624  return cls.fromTable(afwTable)
625 
626  @classmethod
627  def readLsstDefectsFile(cls, filename):
628  """Read defects information from a legacy LSST format text file.
629 
630  Parameters
631  ----------
632  filename : `str`
633  Name of text file containing the defect information.
634 
635  Returns
636  -------
637  defects : `Defects`
638  The defects.
639 
640  Notes
641  -----
642  These defect text files are used as the human readable definitions
643  of defects in calibration data definition repositories. The format
644  is to use four columns defined as follows:
645 
646  x0 : `int`
647  X coordinate of bottom left corner of box.
648  y0 : `int`
649  Y coordinate of bottom left corner of box.
650  width : `int`
651  X extent of the box.
652  height : `int`
653  Y extent of the box.
654 
655  Files of this format were used historically to represent defects
656  in simple text form. Use `Defects.readText` and `Defects.writeText`
657  to use the more modern format.
658  """
659  # Use loadtxt so that ValueError is thrown if the file contains a
660  # non-integer value. genfromtxt converts bad values to -1.
661  defect_array = np.loadtxt(filename,
662  dtype=[("x0", "int"), ("y0", "int"),
663  ("x_extent", "int"), ("y_extent", "int")])
664 
665  return cls(lsst.geom.Box2I(lsst.geom.Point2I(row["x0"], row["y0"]),
666  lsst.geom.Extent2I(row["x_extent"], row["y_extent"]))
667  for row in defect_array)
668 
669  @classmethod
670  def fromFootprintList(cls, fpList):
671  """Compute a defect list from a footprint list, optionally growing
672  the footprints.
673 
674  Parameters
675  ----------
676  fpList : `list` of `lsst.afw.detection.Footprint`
677  Footprint list to process.
678 
679  Returns
680  -------
681  defects : `Defects`
682  List of defects.
683  """
684  return cls(itertools.chain.from_iterable(lsst.afw.detection.footprintToBBoxList(fp)
685  for fp in fpList))
686 
687  @classmethod
688  def fromMask(cls, maskedImage, maskName):
689  """Compute a defect list from a specified mask plane.
690 
691  Parameters
692  ----------
693  maskedImage : `lsst.afw.image.MaskedImage`
694  Image to process.
695  maskName : `str` or `list`
696  Mask plane name, or list of names to convert.
697 
698  Returns
699  -------
700  defects : `Defects`
701  Defect list constructed from masked pixels.
702  """
703  mask = maskedImage.getMask()
704  thresh = lsst.afw.detection.Threshold(mask.getPlaneBitMask(maskName),
705  lsst.afw.detection.Threshold.BITMASK)
706  fpList = lsst.afw.detection.FootprintSet(mask, thresh).getFootprints()
707  return cls.fromFootprintList(fpList)
Defines the fields and offsets for a table.
Definition: Schema.h:50
Encapsulate information about a bad portion of a detector.
Definition: Defect.h:41
def maskPixels(self, maskedImage, maskName="BAD")
Definition: defects.py:252
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
A compact representation of a collection of pixels.
Definition: SpanSet.h:77
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:43
def setMetadata(self, metadata=None)
Definition: defects.py:200
def __setitem__(self, index, value)
Definition: defects.py:151
def __init__(self, defectList=None, metadata=None)
Definition: defects.py:98
def policyToBadRegionList(policyFile)
Definition: defects.py:57
def readLsstDefectsFile(cls, filename)
Definition: defects.py:627
std::vector< lsst::geom::Box2I > footprintToBBoxList(Footprint const &footprint)
Return a list of BBoxs, whose union contains exactly the pixels in the footprint, neither more nor le...
Definition: Footprint.cc:358
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
def insert(self, index, value)
Definition: defects.py:185
static Box2I makeCenteredBox(Point2D const &center, Extent const &size)
Create a box centered as closely as possible on a particular point.
Definition: Box.cc:100
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
An integer coordinate rectangle.
Definition: Box.h:54
def fromMask(cls, maskedImage, maskName)
Definition: defects.py:688
bool strip
Definition: fits.cc:883