LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
registerImage.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008-2013 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 
23 """
24 This module contains a Task to register (align) multiple images.
25 """
26 __all__ = ["RegisterTask", "RegisterConfig"]
27 
28 import math
29 import numpy
30 
31 from lsst.pex.config import Config, Field, ConfigField
32 from lsst.pipe.base import Task, Struct
33 from lsst.meas.astrom.sip import makeCreateWcsWithSip
34 from lsst.afw.math import Warper
35 
36 import lsst.geom as geom
37 import lsst.afw.table as afwTable
38 
39 
41  """Configuration for RegisterTask"""
42  matchRadius = Field(dtype=float, default=1.0, doc="Matching radius (arcsec)", check=lambda x: x > 0)
43  sipOrder = Field(dtype=int, default=4, doc="Order for SIP WCS", check=lambda x: x > 1)
44  sipIter = Field(dtype=int, default=3, doc="Rejection iterations for SIP WCS", check=lambda x: x > 0)
45  sipRej = Field(dtype=float, default=3.0, doc="Rejection threshold for SIP WCS", check=lambda x: x > 0)
46  warper = ConfigField(dtype=Warper.ConfigClass, doc="Configuration for warping")
47 
48 
50  """
51  Task to register (align) multiple images.
52 
53  The 'run' method provides a revised Wcs from matches and fitting sources.
54  Additional methods are provided as a convenience to warp an exposure
55  ('warpExposure') and sources ('warpSources') with the new Wcs.
56  """
57  ConfigClass = RegisterConfig
58 
59  def run(self, inputSources, inputWcs, inputBBox, templateSources):
60  """Register (align) an input exposure to the template
61 
62  The sources must have RA,Dec set, and accurate to within the
63  'matchRadius' of the configuration in order to facilitate source
64  matching. We fit a new Wcs, but do NOT set it in the input exposure.
65 
66  @param inputSources: Sources from input exposure
67  @param inputWcs: Wcs of input exposure
68  @param inputBBox: Bounding box of input exposure
69  @param templateSources: Sources from template exposure
70  @return Struct(matches: Matches between sources,
71  wcs: Wcs for input in frame of template,
72  )
73  """
74  matches = self.matchSourcesmatchSources(inputSources, templateSources)
75  wcs = self.fitWcsfitWcs(matches, inputWcs, inputBBox)
76  return Struct(matches=matches, wcs=wcs)
77 
78  def matchSources(self, inputSources, templateSources):
79  """Match sources between the input and template
80 
81  The order of the input arguments matters (because the later Wcs
82  fitting assumes a particular order).
83 
84  @param inputSources: Source catalog of the input frame
85  @param templateSources: Source of the target frame
86  @return Match list
87  """
88  matches = afwTable.matchRaDec(templateSources, inputSources,
89  self.configconfig.matchRadius*geom.arcseconds)
90  self.loglog.info("Matching within %.1f arcsec: %d matches" % (self.configconfig.matchRadius, len(matches)))
91  self.metadatametadata.set("MATCH_NUM", len(matches))
92  if len(matches) == 0:
93  raise RuntimeError("Unable to match source catalogs")
94  return matches
95 
96  def fitWcs(self, matches, inputWcs, inputBBox):
97  """Fit Wcs to matches
98 
99  The fitting includes iterative sigma-clipping.
100 
101  @param matches: List of matches (first is target, second is input)
102  @param inputWcs: Original input Wcs
103  @param inputBBox: Bounding box of input image
104  @return Wcs
105  """
106  copyMatches = type(matches)(matches)
107  refCoordKey = copyMatches[0].first.getTable().getCoordKey()
108  inCentroidKey = copyMatches[0].second.getTable().getCentroidSlot().getMeasKey()
109  for i in range(self.configconfig.sipIter):
110  sipFit = makeCreateWcsWithSip(copyMatches, inputWcs, self.configconfig.sipOrder, inputBBox)
111  self.loglog.debug("Registration WCS RMS iteration %d: %f pixels",
112  i, sipFit.getScatterInPixels())
113  wcs = sipFit.getNewWcs()
114  dr = [m.first.get(refCoordKey).separation(
115  wcs.pixelToSky(m.second.get(inCentroidKey))).asArcseconds() for
116  m in copyMatches]
117  dr = numpy.array(dr)
118  rms = math.sqrt((dr*dr).mean()) # RMS from zero
119  rms = max(rms, 1.0e-9) # Don't believe any RMS smaller than this
120  self.loglog.debug("Registration iteration %d: rms=%f", i, rms)
121  good = numpy.where(dr < self.configconfig.sipRej*rms)[0]
122  numBad = len(copyMatches) - len(good)
123  self.loglog.debug("Registration iteration %d: rejected %d", i, numBad)
124  if numBad == 0:
125  break
126  copyMatches = type(matches)(copyMatches[i] for i in good)
127 
128  sipFit = makeCreateWcsWithSip(copyMatches, inputWcs, self.configconfig.sipOrder, inputBBox)
129  self.loglog.info("Registration WCS: final WCS RMS=%f pixels from %d matches" %
130  (sipFit.getScatterInPixels(), len(copyMatches)))
131  self.metadatametadata.set("SIP_RMS", sipFit.getScatterInPixels())
132  self.metadatametadata.set("SIP_GOOD", len(copyMatches))
133  self.metadatametadata.set("SIP_REJECTED", len(matches) - len(copyMatches))
134  wcs = sipFit.getNewWcs()
135  return wcs
136 
137  def warpExposure(self, inputExp, newWcs, templateWcs, templateBBox):
138  """Warp input exposure to template frame
139 
140  There are a variety of data attached to the exposure (e.g., PSF, PhotoCalib
141  and other metadata), but we do not attempt to warp these to the template
142  frame.
143 
144  @param inputExp: Input exposure, to be warped
145  @param newWcs: Revised Wcs for input exposure
146  @param templateWcs: Target Wcs
147  @param templateBBox: Target bounding box
148  @return Warped exposure
149  """
150  warper = Warper.fromConfig(self.configconfig.warper)
151  copyExp = inputExp.Factory(inputExp.getMaskedImage(), newWcs)
152  alignedExp = warper.warpExposure(templateWcs, copyExp, destBBox=templateBBox)
153  return alignedExp
154 
155  def warpSources(self, inputSources, newWcs, templateWcs, templateBBox):
156  """Warp sources to the new frame
157 
158  It would be difficult to transform all possible quantities of potential
159  interest between the two frames. We therefore update only the sky and
160  pixel coordinates.
161 
162  @param inputSources: Sources on input exposure, to be warped
163  @param newWcs: Revised Wcs for input exposure
164  @param templateWcs: Target Wcs
165  @param templateBBox: Target bounding box
166  @return Warped sources
167  """
168  alignedSources = inputSources.copy(True)
169  if not isinstance(templateBBox, geom.Box2D):
170  # There is no method Box2I::contains(Point2D)
171  templateBBox = geom.Box2D(templateBBox)
172  table = alignedSources.getTable()
173  coordKey = table.getCoordKey()
174  centroidKey = table.getCentroidSlot().getMeasKey()
175  deleteList = []
176  for i, s in enumerate(alignedSources):
177  oldCentroid = s.get(centroidKey)
178  newCoord = newWcs.pixelToSky(oldCentroid)
179  newCentroid = templateWcs.skyToPixel(newCoord)
180  if not templateBBox.contains(newCentroid):
181  deleteList.append(i)
182  continue
183  s.set(coordKey, newCoord)
184  s.set(centroidKey, newCentroid)
185 
186  for i in reversed(deleteList): # Delete from back so we don't change indices
187  del alignedSources[i]
188 
189  return alignedSources
int max
table::Key< int > type
Definition: Detector.cc:163
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
def matchSources(self, inputSources, templateSources)
def fitWcs(self, matches, inputWcs, inputBBox)
def warpSources(self, inputSources, newWcs, templateWcs, templateBBox)
def run(self, inputSources, inputWcs, inputBBox, templateSources)
def warpExposure(self, inputExp, newWcs, templateWcs, templateBBox)
daf::base::PropertySet * set
Definition: fits.cc:912
std::vector< Match< typename Cat1::Record, typename Cat2::Record > > matchRaDec(Cat1 const &cat1, Cat2 const &cat2, lsst::geom::Angle radius, MatchControl const &mc=MatchControl())
Compute all tuples (s1,s2,d) where s1 belings to cat1, s2 belongs to cat2 and d, the distance between...
Definition: Match.cc:158
CreateWcsWithSip< MatchT > makeCreateWcsWithSip(std::vector< MatchT > const &matches, afw::geom::SkyWcs const &linearWcs, int const order, geom::Box2I const &bbox=geom::Box2I(), int const ngrid=0)
Factory function for CreateWcsWithSip.