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
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 /*****************************************************************************/
28 #include <cassert>
29 #include <string>
30 #include <typeinfo>
31 #include <algorithm>
32 #include "boost/format.hpp"
33 #include "lsst/pex/logging/Trace.h"
34 #include "lsst/pex/exceptions.h"
41 #include "lsst/afw/detection/FootprintArray.cc"
45 
46 namespace lsst {
47 namespace afw {
48 namespace detection {
49 namespace {
50  template<typename T>
51  struct setPixel {
52  setPixel(T val) : _val(val) {}
53 
54  T operator()(T) const {
55  return _val;
56  }
57  private:
58  T _val;
59  };
60 
61  template<>
62  struct setPixel<boost::uint16_t> {
63  typedef boost::uint16_t T;
64 
65  setPixel(T val) : _mask(~val) {}
66 
67  T operator()(T pix) const {
68  pix &= _mask;
69  return pix;
70  }
71  private:
72  T _mask;
73  };
74 }
75 
76 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
78  Footprint const& foot,
80  HeavyFootprintCtrl const *ctrl
81  ) : Footprint(foot),
82  _image(ndarray::allocate(ndarray::makeVector(foot.getNpix()))),
83  _mask(ndarray::allocate(ndarray::makeVector(foot.getNpix()))),
84  _variance(ndarray::allocate(ndarray::makeVector(foot.getNpix())))
85 {
86  HeavyFootprintCtrl ctrl_s = HeavyFootprintCtrl();
87 
88  if (!ctrl) {
89  ctrl = &ctrl_s;
90  }
91 
92  switch (ctrl->getModifySource()) {
94  flattenArray(*this, mimage.getImage()->getArray(), _image, mimage.getXY0());
95  flattenArray(*this, mimage.getMask()->getArray(), _mask, mimage.getXY0());
96  flattenArray(*this, mimage.getVariance()->getArray(), _variance, mimage.getXY0());
97  break;
99  {
100  ImagePixelT const ival = ctrl->getImageVal();
101  MaskPixelT const mval = ctrl->getMaskVal();
102  VariancePixelT const vval = ctrl->getVarianceVal();
103 
104  flattenArray(*this, mimage.getImage()->getArray(), _image,
105  setPixel<ImagePixelT>(ival), mimage.getXY0());
106  flattenArray(*this, mimage.getMask()->getArray(), _mask,
107  setPixel<MaskPixelT>(mval), mimage.getXY0());
108  flattenArray(*this, mimage.getVariance()->getArray(), _variance,
109  setPixel<VariancePixelT>(vval), mimage.getXY0());
110  break;
111  }
112  }
113 }
114 
115 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
117  Footprint const& foot,
118  HeavyFootprintCtrl const* ctrl)
119  : Footprint(foot),
120  _image (ndarray::allocate(ndarray::makeVector(foot.getNpix()))),
121  _mask (ndarray::allocate(ndarray::makeVector(foot.getNpix()))),
122  _variance(ndarray::allocate(ndarray::makeVector(foot.getNpix())))
123 {
124 }
125 
126 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
129  ) const
130 {
131  expandArray(*this, _image, mimage.getImage()->getArray(), mimage.getXY0());
132  expandArray(*this, _mask, mimage.getMask()->getArray(), mimage.getXY0());
133  expandArray(*this, _variance, mimage.getVariance()->getArray(), mimage.getXY0());
134 }
135 
136 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
139 {
140  expandArray(*this, _image, image.getArray(), image.getXY0());
141 }
142 
143 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
144 PTR(HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT>)
145 mergeHeavyFootprints(HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT> const& h1,
146  HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT> const& h2)
147 {
148  // Merge the Footprints (by merging the Spans)
149  PTR(Footprint) foot = mergeFootprints(h1, h2);
150 
151  // Find the union bounding-box
152  geom::Box2I bbox(h1.getBBox());
153  bbox.include(h2.getBBox());
154 
155  // Create union-bb-sized images and insert the heavies
156  image::MaskedImage<ImagePixelT,MaskPixelT,VariancePixelT> im1(bbox);
157  image::MaskedImage<ImagePixelT,MaskPixelT,VariancePixelT> im2(bbox);
158  h1.insert(im1);
159  h2.insert(im2);
160  // Add the pixels
161  im1 += im2;
162 
163  // Build new HeavyFootprint from the merged spans and summed pixels.
164  return PTR(HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT>)
165  (new HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT>(*foot, im1));
166 }
167 
168 
169 template<typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
170 double
171 HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT>::dot(
172  HeavyFootprint<ImagePixelT, MaskPixelT, VariancePixelT> const& rhs
173  ) const
174 {
175  // Require that footprints are sorted in ascending y
176  assert(isNormalized());
177  assert(rhs.isNormalized());
178 
179  // Coordinated cycling through the iterators while juggling the offsets into the arrays
180  typedef typename ndarray::Array<ImagePixelT const, 1, 1>::Iterator ArrayIter;
181  ArrayIter lhsArray = getImageArray().begin(), rhsArray = rhs.getImageArray().begin();
182  SpanList::const_iterator lhsIter = getSpans().begin(), rhsIter = rhs.getSpans().begin();
183  SpanList::const_iterator const lhsEnd = getSpans().end(), rhsEnd = rhs.getSpans().end();
184  double sum = 0.0;
185  while (lhsIter != lhsEnd && rhsIter != rhsEnd) {
186  Span const& lhsSpan = **lhsIter, rhsSpan = **rhsIter;
187  int const yLhs = lhsSpan.getY(), yRhs = rhsSpan.getY();
188  if (yLhs == yRhs) {
189  int const x0Lhs = lhsSpan.getX0(), x1Lhs = lhsSpan.getX1();
190  int const x0Rhs = rhsSpan.getX0(), x1Rhs = rhsSpan.getX1();
191  int const xMin = std::max(x0Lhs, x0Rhs), xMax = std::min(x1Lhs, x1Rhs);
192  if (xMin <= xMax) {
193  lhsArray += xMin - x0Lhs;
194  rhsArray += xMin - x0Rhs;
195  for (int x = xMin; x <= xMax; ++x, ++lhsArray, ++rhsArray) {
196  sum += (*lhsArray)*(*rhsArray);
197  }
198  // Rewind to the start of the span, for easier sync between spans and arrays
199  lhsArray -= xMax + 1 - x0Lhs;
200  rhsArray -= xMax + 1 - x0Rhs;
201  }
202  if (x1Lhs <= x1Rhs) {
203  lhsArray += lhsSpan.getWidth();
204  ++lhsIter;
205  } else {
206  rhsArray += rhsSpan.getWidth();
207  ++rhsIter;
208  }
209  continue;
210  } else if (yLhs < yRhs) {
211  while (lhsIter != lhsEnd && (*lhsIter)->getY() < yRhs) {
212  lhsArray += (*lhsIter)->getWidth();
213  ++lhsIter;
214  }
215  continue;
216  } else { // yLhs > yRhs
217  while (rhsIter != rhsEnd && (*rhsIter)->getY() < yLhs) {
218  rhsArray += (*rhsIter)->getWidth();
219  ++rhsIter;
220  }
221  continue;
222  }
223  }
224  return sum;
225 }
226 
227 
228 /************************************************************************************************************/
229 // Persistence (using afw::table::io)
230 //
231 
232 namespace {
233 
234 // Schema and Keys used to persist the pixels of a HeavyFootprint (Spans and Peaks are handled by the
235 // Footprint base class). This is a singleton, but a different one for each template instantiation.
236 template <typename ImagePixelT,
237  typename MaskPixelT=image::MaskPixel,
238  typename VariancePixelT=image::VariancePixel>
239 struct HeavyFootprintPersistenceHelper {
240  afw::table::Schema schema;
241  afw::table::Key< afw::table::Array<ImagePixelT> > image;
242  afw::table::Key< afw::table::Array<MaskPixelT> > mask;
243  afw::table::Key< afw::table::Array<VariancePixelT> > variance;
244 
245  static HeavyFootprintPersistenceHelper const & get() {
246  static HeavyFootprintPersistenceHelper const instance;
247  return instance;
248  }
249 
250 private:
251 
252  HeavyFootprintPersistenceHelper() :
253  schema(),
254  image(schema.addField< afw::table::Array<ImagePixelT> >(
255  "image", "image pixels for HeavyFootprint", "dn"
256  )),
257  mask(schema.addField< afw::table::Array<MaskPixelT> >(
258  "mask", "mask pixels for HeavyFootprint"
259  )),
260  variance(schema.addField< afw::table::Array<VariancePixelT> >(
261  "variance", "variance pixels for HeavyFootprint", "dn^2"
262  ))
263  {
264  schema.getCitizen().markPersistent();
265  }
266 
267 };
268 
269 // These suffix-computing structs are used to compute the string name associated with a HeavyFootprint
270 // for Persistence.
271 // We don't instantiate HeavyFootprints with anything other than defaults for Mask and Variance, so we
272 // don't bother figuring out what suffixes to use for them for now. If we change that, we just need
273 // to add more explicit specializations of this template.
274 template <typename ImagePixelT,
275  typename MaskPixelT=image::MaskPixel,
276  typename VariancePixelT=image::VariancePixel>
277 struct ComputeSuffix;
278 template <> struct ComputeSuffix<boost::uint16_t> { static std::string apply() { return "U"; } };
279 template <> struct ComputeSuffix<float> { static std::string apply() { return "F"; } };
280 template <> struct ComputeSuffix<double> { static std::string apply() { return "D"; } };
281 template <> struct ComputeSuffix<int> { static std::string apply() { return "I"; } };
282 
283 } // anonymous
284 
285 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
287  return "HeavyFootprint" + ComputeSuffix<ImagePixelT,MaskPixelT,VariancePixelT>::apply();
288 }
289 
290 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
291 void HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT>::write(OutputArchiveHandle & handle) const {
292  HeavyFootprintPersistenceHelper<ImagePixelT,MaskPixelT,VariancePixelT> const & keys =
293  HeavyFootprintPersistenceHelper<ImagePixelT,MaskPixelT,VariancePixelT>::get();
294  // delegate to Footprint::write to handle spans and peaks
295  Footprint::write(handle);
296  // add one more catalog for pixel values
297  afw::table::BaseCatalog cat = handle.makeCatalog(keys.schema);
298  PTR(afw::table::BaseRecord) record = cat.addNew();
299  // We could deep-copy the arrays instead of const-casting them, which might be marginally safer,
300  // but we always save an OutputArchive to disk immediately after we create it, so there's really
301  // no chance that we could get the HeavyFootprint in trouble by having this view modified.
302  record->set(keys.image, ndarray::const_array_cast<ImagePixelT>(getImageArray()));
303  record->set(keys.mask, ndarray::const_array_cast<MaskPixelT>(getMaskArray()));
304  record->set(keys.variance, ndarray::const_array_cast<VariancePixelT>(getVarianceArray()));
305  handle.saveCatalog(cat);
306 }
307 
308 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
309 class HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT>::Factory :
310  public afw::table::io::PersistableFactory
311 {
312 public:
313 
314  explicit Factory(std::string const & name) : afw::table::io::PersistableFactory(name) {}
315 
316  virtual PTR(afw::table::io::Persistable)
317  read(InputArchive const & archive, CatalogVector const & catalogs) const {
318  HeavyFootprintPersistenceHelper<ImagePixelT,MaskPixelT,VariancePixelT> const & keys =
319  HeavyFootprintPersistenceHelper<ImagePixelT,MaskPixelT,VariancePixelT>::get();
320  LSST_ARCHIVE_ASSERT(catalogs.size() == 3u);
321  PTR(HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT>) result(
322  new HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT>()
323  );
324  result->readSpans(catalogs[0]); // these read methods are inherited from Footprint
325  result->readPeaks(catalogs[1]);
326  afw::table::BaseRecord const & record = catalogs[2].front();
327  result->_image = ndarray::const_array_cast<ImagePixelT>(record.get(keys.image));
328  result->_mask = ndarray::const_array_cast<MaskPixelT>(record.get(keys.mask));
329  result->_variance = ndarray::const_array_cast<VariancePixelT>(record.get(keys.variance));
330  return result;
331  }
332 
333  static Factory registration;
334 
335 };
336 
337 // initialize static instance, registering the factory with the persistence mechanism at the same time
338 template <typename ImagePixelT, typename MaskPixelT, typename VariancePixelT>
339 typename HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT>::Factory
340 HeavyFootprint<ImagePixelT,MaskPixelT,VariancePixelT>::Factory::registration(
341  "HeavyFootprint" + ComputeSuffix<ImagePixelT,MaskPixelT,VariancePixelT>::apply()
342 );
343 
344 /************************************************************************************************************/
345 //
346 // Explicit instantiations
347 // \cond
348 //
349 //
350 #define INSTANTIATE(TYPE) \
351  template class HeavyFootprint<TYPE>; \
352  template PTR(HeavyFootprint<TYPE>) mergeHeavyFootprints<TYPE>( \
353  HeavyFootprint<TYPE> const&, HeavyFootprint<TYPE> const&);
354 
355 INSTANTIATE(boost::uint16_t);
356 INSTANTIATE(double);
357 INSTANTIATE(float);
358 INSTANTIATE(int);
359 
360 }}}
361 // \endcond
void flattenArray(Footprint const &fp, ndarray::Array< T, N, C > const &src, ndarray::Array< U, N-1, D > const &dest, lsst::afw::geom::Point2I const &xy0=lsst::afw::geom::Point2I())
Flatten the first two dimensions of an array.
boost::uint16_t MaskPixel
ImagePtr getImage(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:869
void insert(lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage) const
#define PTR(...)
Definition: base.h:41
Represent a set of pixels of an arbitrary shape and size.
table::Key< std::string > name
Definition: ApCorrMap.cc:71
void expandArray(Footprint const &fp, ndarray::Array< T, N, C > const &src, ndarray::Array< U, N+1, D > const &dest, lsst::afw::geom::Point2I const &xy0=lsst::afw::geom::Point2I())
expand the first dimension of an array
CatalogT< BaseRecord > BaseCatalog
Definition: fwd.h:61
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:47
Array< T_, N, C > const_array_cast(Array< T, N, C > const &array)
Definition: casts.h:76
lsst::afw::detection::Footprint Footprint
Definition: Source.h:61
definition of the Trace messaging facilities
Represent a set of pixels of an arbitrary shape and size, including values for those pixels; a HeavyF...
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
Definition: operators.h:1250
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
VariancePtr getVariance(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s variance.
Definition: MaskedImage.h:890
virtual void write(OutputArchiveHandle &handle) const
Write the object to one or more catalogs.
virtual void write(OutputArchiveHandle &handle) const
Write the object to one or more catalogs.
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:77
Vector< T, N > makeVector(T v1, T v2,..., T vN)
Variadic constructor for Vector.
float VariancePixel
! default type for Masks and MaskedImage Masks
geom::Point2I getXY0() const
Definition: MaskedImage.h:929
tbl::Schema schema
int x
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
#define INSTANTIATE(T)
Control Footprint-related algorithms.
MaskPtr getMask(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s mask.
Definition: MaskedImage.h:879
virtual std::string getPersistenceName() const
Return the unique name used to persist this object and look up its factory.
geom::Point2I getXY0() const
Definition: Image.h:264
Implementation of the Class MaskedImage.
ExpressionTraits< Derived >::Iterator Iterator
Nested expression or element iterator.
boost::shared_ptr< Footprint > mergeFootprints(Footprint const &foot1, Footprint const &foot2)
ImageT val
Definition: CR.cc:154
boost::shared_ptr< HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > > mergeHeavyFootprints(HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > const &h1, HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > const &h2)
Include files required for standard LSST Exception handling.