LSSTApplications  18.0.0+106,18.0.0+50,19.0.0,19.0.0+1,19.0.0+10,19.0.0+11,19.0.0+13,19.0.0+17,19.0.0+2,19.0.0-1-g20d9b18+6,19.0.0-1-g425ff20,19.0.0-1-g5549ca4,19.0.0-1-g580fafe+6,19.0.0-1-g6fe20d0+1,19.0.0-1-g7011481+9,19.0.0-1-g8c57eb9+6,19.0.0-1-gb5175dc+11,19.0.0-1-gdc0e4a7+9,19.0.0-1-ge272bc4+6,19.0.0-1-ge3aa853,19.0.0-10-g448f008b,19.0.0-12-g6990b2c,19.0.0-2-g0d9f9cd+11,19.0.0-2-g3d9e4fb2+11,19.0.0-2-g5037de4,19.0.0-2-gb96a1c4+3,19.0.0-2-gd955cfd+15,19.0.0-3-g2d13df8,19.0.0-3-g6f3c7dc,19.0.0-4-g725f80e+11,19.0.0-4-ga671dab3b+1,19.0.0-4-gad373c5+3,19.0.0-5-ga2acb9c+2,19.0.0-5-gfe96e6c+2,w.2020.01
LSSTDataManagementBasePackage
HeavyFootprint.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2008, 2009, 2010 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 #include <cstdint>
24 #include <cassert>
25 #include <string>
26 #include <typeinfo>
27 #include <algorithm>
28 #include "boost/format.hpp"
29 #include "lsst/pex/exceptions.h"
40 
41 namespace lsst {
42 namespace afw {
43 namespace detection {
44 namespace {
45 
46 template <typename T>
47 struct FlattenWithSetter {
48  FlattenWithSetter(T val) : _val(val) {}
49 
50  void operator()(lsst::geom::Point2I const& point, T& out, T& in) {
51  out = in;
52  in = _val;
53  }
54 
55 private:
56  T _val;
57 };
58 
59 template <>
60 struct FlattenWithSetter<lsst::afw::image::MaskPixel> {
62  FlattenWithSetter(T val) : _mask(~val) {}
63 
64  void operator()(lsst::geom::Point2I const& point, T& out, T& in) {
65  out = in;
66  in &= _mask;
67  }
68 
69 private:
70  T _mask;
71 };
72 } // namespace
73 
74 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
77  HeavyFootprintCtrl const* ctrl)
78  : Footprint(foot),
79  _image(ndarray::allocate(ndarray::makeVector(foot.getArea()))),
80  _mask(ndarray::allocate(ndarray::makeVector(foot.getArea()))),
81  _variance(ndarray::allocate(ndarray::makeVector(foot.getArea()))) {
83 
84  if (!ctrl) {
85  ctrl = &ctrl_s;
86  }
87 
88  switch (ctrl->getModifySource()) {
90  getSpans()->flatten(_image, mimage.getImage()->getArray(), mimage.getXY0());
91  getSpans()->flatten(_mask, mimage.getMask()->getArray(), mimage.getXY0());
92  getSpans()->flatten(_variance, mimage.getVariance()->getArray(), mimage.getXY0());
93  break;
95  ImagePixelT const ival = ctrl->getImageVal();
96  MaskPixelT const mval = ctrl->getMaskVal();
97  VariancePixelT const vval = ctrl->getVarianceVal();
98 
99  getSpans()->applyFunctor(FlattenWithSetter<ImagePixelT>(ival), ndarray::ndFlat(_image),
100  ndarray::ndImage(mimage.getImage()->getArray(), mimage.getXY0()));
101  getSpans()->applyFunctor(FlattenWithSetter<MaskPixelT>(mval), ndarray::ndFlat(_mask),
102  ndarray::ndImage(mimage.getMask()->getArray(), mimage.getXY0()));
103  getSpans()->applyFunctor(FlattenWithSetter<VariancePixelT>(vval), ndarray::ndFlat(_variance),
104  ndarray::ndImage(mimage.getVariance()->getArray(), mimage.getXY0()));
105  break;
106  }
107  }
108 }
109 
110 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
112  HeavyFootprintCtrl const* ctrl)
113  : Footprint(foot),
114  _image(ndarray::allocate(ndarray::makeVector(foot.getArea()))),
115  _mask(ndarray::allocate(ndarray::makeVector(foot.getArea()))),
116  _variance(ndarray::allocate(ndarray::makeVector(foot.getArea()))) {}
117 
118 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
121  getSpans()->unflatten(mimage.getImage()->getArray(), _image, mimage.getXY0());
122  getSpans()->unflatten(mimage.getMask()->getArray(), _mask, mimage.getXY0());
123  getSpans()->unflatten(mimage.getVariance()->getArray(), _variance, mimage.getXY0());
124 }
125 
126 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
128  getSpans()->unflatten(image.getArray(), _image, image.getXY0());
129 }
130 
131 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
135  // Merge the Footprints (by merging the Spans)
137 
138  // Find the union bounding-box
140  bbox.include(h2.getBBox());
141 
142  // Create union-bb-sized images and insert the heavies
145  h1.insert(im1);
146  h2.insert(im2);
147  // Add the pixels
148  im1 += im2;
149 
150  // Build new HeavyFootprint from the merged spans and summed pixels.
151  return std::make_shared<HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT>>(*foot, im1);
152 }
153 
154 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
157  // Coordinated cycling through the iterators while juggling the offsets into the arrays
158  typedef typename ndarray::Array<ImagePixelT const, 1, 1>::Iterator ArrayIter;
159  ArrayIter lhsArray = getImageArray().begin(), rhsArray = rhs.getImageArray().begin();
160  auto lhsIter = getSpans()->begin(), rhsIter = rhs.getSpans()->begin();
161  auto const lhsEnd = getSpans()->end(), rhsEnd = rhs.getSpans()->end();
162  double sum = 0.0;
163  while (lhsIter != lhsEnd && rhsIter != rhsEnd) {
164  geom::Span const &lhsSpan = *lhsIter, rhsSpan = *rhsIter;
165  int const yLhs = lhsSpan.getY(), yRhs = rhsSpan.getY();
166  if (yLhs == yRhs) {
167  int const x0Lhs = lhsSpan.getX0(), x1Lhs = lhsSpan.getX1();
168  int const x0Rhs = rhsSpan.getX0(), x1Rhs = rhsSpan.getX1();
169  int const xMin = std::max(x0Lhs, x0Rhs), xMax = std::min(x1Lhs, x1Rhs);
170  if (xMin <= xMax) {
171  lhsArray += xMin - x0Lhs;
172  rhsArray += xMin - x0Rhs;
173  for (int x = xMin; x <= xMax; ++x, ++lhsArray, ++rhsArray) {
174  sum += (*lhsArray) * (*rhsArray);
175  }
176  // Rewind to the start of the span, for easier sync between spans and arrays
177  lhsArray -= xMax + 1 - x0Lhs;
178  rhsArray -= xMax + 1 - x0Rhs;
179  }
180  if (x1Lhs <= x1Rhs) {
181  lhsArray += lhsSpan.getWidth();
182  ++lhsIter;
183  } else {
184  rhsArray += rhsSpan.getWidth();
185  ++rhsIter;
186  }
187  continue;
188  } else if (yLhs < yRhs) {
189  while (lhsIter != lhsEnd && lhsIter->getY() < yRhs) {
190  lhsArray += lhsIter->getWidth();
191  ++lhsIter;
192  }
193  continue;
194  } else { // yLhs > yRhs
195  while (rhsIter != rhsEnd && rhsIter->getY() < yLhs) {
196  rhsArray += rhsIter->getWidth();
197  ++rhsIter;
198  }
199  continue;
200  }
201  }
202  return sum;
203 }
204 
205 // Persistence (using afw::table::io)
206 //
207 
208 namespace {
209 
210 // Schema and Keys used to persist the pixels of a HeavyFootprint (Spans and Peaks are handled by the
211 // Footprint base class). This is a singleton, but a different one for each template instantiation.
212 template <typename ImagePixelT, typename MaskPixelT = image::MaskPixel,
213  typename VariancePixelT = image::VariancePixel>
214 struct HeavyFootprintPersistenceHelper {
215  afw::table::Schema schema;
216  afw::table::Key<afw::table::Array<ImagePixelT>> image;
217  afw::table::Key<afw::table::Array<MaskPixelT>> mask;
218  afw::table::Key<afw::table::Array<VariancePixelT>> variance;
219 
220  static HeavyFootprintPersistenceHelper const& get() {
221  static HeavyFootprintPersistenceHelper const instance;
222  return instance;
223  }
224 
225 private:
226  HeavyFootprintPersistenceHelper()
227  : schema(),
228  image(schema.addField<afw::table::Array<ImagePixelT>>(
229  "image", "image pixels for HeavyFootprint", "count")),
230  mask(schema.addField<afw::table::Array<MaskPixelT>>("mask", "mask pixels for HeavyFootprint")),
231  variance(schema.addField<afw::table::Array<VariancePixelT>>(
232  "variance", "variance pixels for HeavyFootprint", "count^2")) {}
233 };
234 
235 // These suffix-computing structs are used to compute the string name associated with a HeavyFootprint
236 // for Persistence.
237 // We don't instantiate HeavyFootprints with anything other than defaults for Mask and Variance, so we
238 // don't bother figuring out what suffixes to use for them for now. If we change that, we just need
239 // to add more explicit specializations of this template.
240 template <typename ImagePixelT, typename MaskPixelT = image::MaskPixel,
241  typename VariancePixelT = image::VariancePixel>
242 struct ComputeSuffix;
243 template <>
244 struct ComputeSuffix<std::uint16_t> {
245  static std::string apply() { return "U"; }
246 };
247 template <>
248 struct ComputeSuffix<float> {
249  static std::string apply() { return "F"; }
250 };
251 template <>
252 struct ComputeSuffix<double> {
253  static std::string apply() { return "D"; }
254 };
255 template <>
256 struct ComputeSuffix<int> {
257  static std::string apply() { return "I"; }
258 };
259 
260 } // namespace
261 
262 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
264  return "HeavyFootprint" + ComputeSuffix<ImagePixelT, MaskPixelT, VariancePixelT>::apply();
265 }
266 
267 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
269  HeavyFootprintPersistenceHelper<ImagePixelT, MaskPixelT, VariancePixelT> const& keys =
270  HeavyFootprintPersistenceHelper<ImagePixelT, MaskPixelT, VariancePixelT>::get();
271  // delegate to Footprint::write to handle spans and peaks
272  Footprint::write(handle);
273  // add one more catalog for pixel values
274  afw::table::BaseCatalog cat = handle.makeCatalog(keys.schema);
276  // We could deep-copy the arrays instead of const-casting them, which might be marginally safer,
277  // but we always save an OutputArchive to disk immediately after we create it, so there's really
278  // no chance that we could get the HeavyFootprint in trouble by having this view modified.
279  record->set(keys.image, ndarray::const_array_cast<ImagePixelT>(getImageArray()));
280  record->set(keys.mask, ndarray::const_array_cast<MaskPixelT>(getMaskArray()));
281  record->set(keys.variance, ndarray::const_array_cast<VariancePixelT>(getVarianceArray()));
282  handle.saveCatalog(cat);
283 }
284 
285 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
286 class HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT>::Factory
288 public:
289  explicit Factory(std::string const& name) : afw::table::io::PersistableFactory(name) {}
290 
292  CatalogVector const& catalogs) const override {
293  HeavyFootprintPersistenceHelper<ImagePixelT, MaskPixelT, VariancePixelT> const& keys =
294  HeavyFootprintPersistenceHelper<ImagePixelT, MaskPixelT, VariancePixelT>::get();
295  HeavyFootprintPersistenceHelper<ImagePixelT, std::uint16_t, VariancePixelT> const& legacyKeys =
296  HeavyFootprintPersistenceHelper<ImagePixelT, std::uint16_t, VariancePixelT>::get();
297  LSST_ARCHIVE_ASSERT(catalogs.size() == 3u);
298 
299  // Read in the SpanSet into a new Footprint object
300  std::shared_ptr<Footprint> loadedFootprint = readSpanSet(catalogs[0], archive);
301  // Now read in the PeakCatalog records
302  readPeaks(catalogs[1], *loadedFootprint);
303  afw::table::BaseRecord const& record = catalogs[2].front();
304 
305  // Create the HeavyFootprint from the above Footprint
306  auto result =
307  std::make_shared<HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT>>(*loadedFootprint);
308  result->_image = ndarray::const_array_cast<ImagePixelT>(record.get(keys.image));
309 
310  // Handle legacy Masks prior to change to int32
311  if (catalogs[2].getSchema() == legacyKeys.schema) {
312  auto legacyMask = ndarray::const_array_cast<std::uint16_t>(record.get(legacyKeys.mask));
313  result->_mask.deep() = legacyMask;
314  } else {
315  result->_mask = ndarray::const_array_cast<MaskPixelT>(record.get(keys.mask));
316  }
317  result->_variance = ndarray::const_array_cast<VariancePixelT>(record.get(keys.variance));
318  return result;
319  }
320 
322 };
323 
324 // initialize static instance, registering the factory with the persistence mechanism at the same time
325 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
328  "HeavyFootprint" + ComputeSuffix<ImagePixelT, MaskPixelT, VariancePixelT>::apply());
329 } // namespace detection
330 
331 //
332 // Explicit instantiations
333 //
334 //
335 #define INSTANTIATE(TYPE) \
336  template std::shared_ptr<detection::HeavyFootprint<TYPE>> \
337  table::io::PersistableFacade<detection::HeavyFootprint<TYPE>>::dynamicCast( \
338  std::shared_ptr<table::io::Persistable> const&); \
339  template class detection::HeavyFootprint<TYPE>; \
340  template std::shared_ptr<detection::HeavyFootprint<TYPE>> detection::mergeHeavyFootprints<TYPE>( \
341  detection::HeavyFootprint<TYPE> const&, detection::HeavyFootprint<TYPE> const&);
342 
344 INSTANTIATE(double);
345 INSTANTIATE(float);
346 INSTANTIATE(int);
347 } // namespace afw
348 } // namespace lsst
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage&#39;s variance.
Definition: MaskedImage.h:1090
afw::table::Key< afw::table::Array< VariancePixelT > > variance
afw::table::Key< afw::table::Array< ImagePixelT > > image
afw::table::Key< afw::table::Array< MaskPixelT > > mask
An object passed to Persistable::write to allow it to persist itself.
A control object for HeavyFootprints.
Definition: FootprintCtrl.h:99
T front(T... args)
std::shared_ptr< afw::table::io::Persistable > read(InputArchive const &archive, CatalogVector const &catalogs) const override
Construct a new object from the given InputArchive and vector of catalogs.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
void write(OutputArchiveHandle &handle) const override
Write an instance of a Footprint to an output Archive.
A base class for factory classes used to reconstruct objects from records.
Definition: Persistable.h:228
ImageT val
Definition: CR.cc:146
static void readPeaks(afw::table::BaseCatalog const &, Footprint &)
Static method used to unpersist the PeakCatalog member of the Footprint class.
Definition: Footprint.cc:296
lsst::geom::Point2I getXY0() const
Return the image&#39;s origin.
Definition: MaskedImage.h:1121
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
ndarray::Array< MaskPixelT, 1, 1 > getMaskArray()
A set of pixels in an Image, including those pixels&#39; actual values.
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:1057
STL class.
lsst::geom::Box2I getBBox() const
Return the Footprint&#39;s bounding box.
Definition: Footprint.h:208
T min(T... args)
FastFinder::Iterator Iterator
Definition: FastFinder.cc:179
A base class for image defects.
static std::unique_ptr< Footprint > readSpanSet(afw::table::BaseCatalog const &, afw::table::io::InputArchive const &)
Static method used to unpersist the SpanSet member of the Footprint class.
Definition: Footprint.cc:272
afw::table::Schema schema
ndarray::Array< VariancePixelT, 1, 1 > getVarianceArray()
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
details::FlatNdGetter< T, inA, inB > ndFlat(ndarray::Array< T, inA, inB > const &array)
Marks a ndarray to be interpreted as a 1D vector when applying a functor from a SpanSet.
std::size_t getArea() const
Return the number of pixels in this Footprint.
Definition: Footprint.h:173
std::shared_ptr< Footprint > mergeFootprints(Footprint const &footprint1, Footprint const &footprint2)
Merges two Footprints – appends their peaks, and unions their spans, returning a new Footprint...
Definition: Footprint.cc:327
std::string getPersistenceName() const override
Return the name correspoinging ot the persistence type.
std::int32_t MaskPixel
default type for Masks and MaskedImage Masks
T max(T... args)
Class to describe the properties of a detected object from an image.
Definition: Footprint.h:62
double x
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
std::shared_ptr< HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > > mergeHeavyFootprints(HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > const &h1, HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > const &h2)
Sum the two given HeavyFootprints h1 and h2, returning a HeavyFootprint with the union footprint...
MaskPtr getMask() const
Return a (shared_ptr to) the MaskedImage&#39;s mask.
Definition: MaskedImage.h:1069
lsst::geom::Point2I getXY0() const
Return the image&#39;s origin.
Definition: ImageBase.h:360
T size(T... args)
A vector of catalogs used by Persistable.
Definition: CatalogVector.h:29
void write(OutputArchiveHandle &handle) const override
Write an instance of a Footprint to an output Archive.
Definition: Footprint.cc:257
Base class for all records.
Definition: BaseRecord.h:31
void insert(lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage) const
Replace all the pixels in the image with the values in the HeavyFootprint.
ndarray::Array< ImagePixelT, 1, 1 > getImageArray()
std::shared_ptr< geom::SpanSet > getSpans() const
Return a shared pointer to the SpanSet.
Definition: Footprint.h:115
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
double dot(HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > const &other) const
Dot product between HeavyFootprints.
A multi-catalog archive object used to load table::io::Persistable objects.
Definition: InputArchive.h:31
#define INSTANTIATE(TYPE)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
details::ImageNdGetter< T, inA, inB > ndImage(ndarray::Array< T, inA, inB > const &array, lsst::geom::Point2I xy0=lsst::geom::Point2I())
Marks a ndarray to be interpreted as an image when applying a functor from a SpanSet.
float VariancePixel
default type for MaskedImage variance images
An integer coordinate rectangle.
Definition: Box.h:55
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
py::object result
Definition: _schema.cc:429
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:485
HeavyFootprint()=default
Default constructor for HeavyFootprint.