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
Image.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008-2016 AURA/LSST.
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 
25 /*
26  * Implementation for ImageBase and Image
27  */
28 #include <cstdint>
29 #include <iostream>
30 #include <functional>
31 #include <type_traits>
32 #include "boost/mpl/vector.hpp"
33 #pragma clang diagnostic push
34 #pragma clang diagnostic ignored "-Wunused-variable"
35 #pragma clang diagnostic pop
36 #include "boost/format.hpp"
37 #include "boost/filesystem/path.hpp"
38 
39 #include "boost/version.hpp"
40 #if BOOST_VERSION < 106900
41 #include "boost/gil/gil_all.hpp"
42 #else
43 #include "boost/gil.hpp"
44 #endif
45 
46 #include "lsst/pex/exceptions.h"
47 #include "lsst/afw/geom/wcsUtils.h"
48 #include "lsst/afw/image/Image.h"
50 #include "lsst/afw/fits.h"
52 
53 namespace lsst {
54 namespace afw {
55 namespace image {
56 
57 template <typename PixelT>
58 typename ImageBase<PixelT>::_view_t ImageBase<PixelT>::_allocateView(lsst::geom::Extent2I const& dimensions,
59  Manager::Ptr& manager) {
60  if (dimensions.getX() < 0 || dimensions.getY() < 0) {
62  str(boost::format("Both width and height must be non-negative: %d, %d") %
63  dimensions.getX() % dimensions.getY()));
64  }
65  if (dimensions.getX() != 0 && dimensions.getY() > std::numeric_limits<int>::max() / dimensions.getX()) {
67  str(boost::format("Image dimensions (%d x %d) too large; int overflow detected.") %
68  dimensions.getX() % dimensions.getY()));
69  }
71  ndarray::SimpleManager<PixelT>::allocate(dimensions.getX() * dimensions.getY());
72  manager = r.first;
73  return boost::gil::interleaved_view(dimensions.getX(), dimensions.getY(),
74  (typename _view_t::value_type*)r.second,
75  dimensions.getX() * sizeof(PixelT));
76 }
77 template <typename PixelT>
78 typename ImageBase<PixelT>::_view_t ImageBase<PixelT>::_makeSubView(lsst::geom::Extent2I const& dimensions,
79  lsst::geom::Extent2I const& offset,
80  const _view_t& view) {
81  if (offset.getX() < 0 || offset.getY() < 0 || offset.getX() + dimensions.getX() > view.width() ||
82  offset.getY() + dimensions.getY() > view.height()) {
83  throw LSST_EXCEPT(
86  "Box2I(Point2I(%d,%d),lsst::geom::Extent2I(%d,%d)) doesn't fit in image %dx%d") %
87  offset.getX() % offset.getY() % dimensions.getX() % dimensions.getY() % view.width() %
88  view.height())
89  .str());
90  }
91  if (dimensions.getX() == 0 && dimensions.getY() == 0
92  && view.width() == 0 && view.height() == 0) {
93  // Getting a zero-extent subview of a zero-extent view returns itself
94  return view;
95  } else {
96  return boost::gil::subimage_view(view, offset.getX(), offset.getY(), dimensions.getX(),
97  dimensions.getY());
98  }
99 }
100 
101 template <typename PixelT>
103  : _origin(0, 0), _manager(), _gilView(_allocateView(dimensions, _manager)) {}
104 
105 template <typename PixelT>
107  : _origin(bbox.getMin()), _manager(), _gilView(_allocateView(bbox.getDimensions(), _manager)) {}
108 
109 template <typename PixelT>
110 ImageBase<PixelT>::ImageBase(ImageBase const& rhs, bool const deep
111 
112  )
113  : _origin(rhs._origin), _manager(rhs._manager), _gilView(rhs._gilView) {
114  if (deep) {
115  ImageBase tmp(getBBox());
116  tmp.assign(*this); // now copy the pixels
117  swap(tmp);
118  }
119 }
120 // Delegate to copy-constructor for backwards compatibility
121 template <typename PixelT>
123 
124 template <typename PixelT>
126  bool const deep
127 
128  )
129  : _origin((origin == PARENT) ? bbox.getMin() : rhs._origin + lsst::geom::Extent2I(bbox.getMin())),
130  _manager(rhs._manager), // reference counted pointer, don't copy pixels
131  _gilView(_makeSubView(bbox.getDimensions(), _origin - rhs._origin, rhs._gilView)) {
132  if (deep) {
133  ImageBase tmp(getBBox());
134  tmp.assign(*this); // now copy the pixels
135  swap(tmp);
136  }
137 }
138 
139 template <typename PixelT>
140 ImageBase<PixelT>::ImageBase(Array const& array, bool deep, lsst::geom::Point2I const& xy0)
141  : _origin(xy0),
142  _manager(array.getManager()),
143  _gilView(boost::gil::interleaved_view(array.template getSize<1>(), array.template getSize<0>(),
144  (typename _view_t::value_type*)array.getData(),
145  array.template getStride<0>() * sizeof(PixelT))) {
146  if (deep) {
147  ImageBase tmp(*this, true);
148  swap(tmp);
149  }
150 }
151 
152 template <typename PixelT>
154  ImageBase tmp(rhs);
155  swap(tmp); // See Meyers, Effective C++, Item 11
156 
157  return *this;
158 }
159 // Delegate to copy-assignment for backwards compatibility
160 template <typename PixelT>
162  return *this = rhs;
163 }
164 
165 template <typename PixelT>
167  auto lhsDim = bbox.isEmpty() ? getDimensions() : bbox.getDimensions();
168  if (lhsDim != rhs.getDimensions()) {
170  (boost::format("Dimension mismatch: %dx%d v. %dx%d") % lhsDim.getX() %
171  lhsDim.getY() % rhs.getWidth() % rhs.getHeight())
172  .str());
173  }
174  if (bbox.isEmpty()) {
175  copy_pixels(rhs._gilView, _gilView);
176  } else {
177  auto lhsOff = (origin == PARENT) ? bbox.getMin() - _origin : lsst::geom::Extent2I(bbox.getMin());
178  auto lhsGilView = _makeSubView(lhsDim, lhsOff, _gilView);
179  copy_pixels(rhs._gilView, lhsGilView);
180  }
181 }
182 
183 template <typename PixelT>
185  return const_cast<typename ImageBase<PixelT>::PixelReference>(
186  static_cast<typename ImageBase<PixelT>::PixelConstReference>(_gilView(x, y)[0]));
187 }
188 
189 template <typename PixelT>
191  CheckIndices const& check) {
192  if (check && (x < 0 || x >= getWidth() || y < 0 || y >= getHeight())) {
194  (boost::format("Index (%d, %d) is out of range [0--%d], [0--%d]") % x % y %
195  (getWidth() - 1) % (getHeight() - 1))
196  .str());
197  }
198 
199  return const_cast<typename ImageBase<PixelT>::PixelReference>(
200  static_cast<typename ImageBase<PixelT>::PixelConstReference>(_gilView(x, y)[0]));
201 }
202 
203 template <typename PixelT>
205  return _gilView(x, y)[0];
206 }
207 
208 template <typename PixelT>
210  int x, int y, CheckIndices const& check) const {
211  if (check && (x < 0 || x >= getWidth() || y < 0 || y >= getHeight())) {
213  (boost::format("Index (%d, %d) is out of range [0--%d], [0--%d]") % x % y %
214  (this->getWidth() - 1) % (this->getHeight() - 1))
215  .str());
216  }
217 
218  return _gilView(x, y)[0];
219 }
220 
221 template <typename PixelT>
223  ImageOrigin origin) {
224  int x = index.getX();
225  int y = index.getY();
226  if (origin == PARENT) {
227  x -= getX0();
228  y -= getY0();
229  }
230  return _gilView(x, y)[0];
231 }
232 
233 template <typename PixelT>
235  ImageOrigin origin) const {
236  int x = index.getX();
237  int y = index.getY();
238  if (origin == PARENT) {
239  x -= getX0();
240  y -= getY0();
241  }
242  return _gilView(x, y)[0];
243 }
244 
245 template <typename PixelT>
247  using std::swap; // See Meyers, Effective C++, Item 25
248 
249  swap(_manager, rhs._manager); // just swapping the pointers
250  swap(_gilView, rhs._gilView);
251  swap(_origin, rhs._origin);
252 }
253 
254 template <typename PixelT>
256  a.swap(b);
257 }
258 
259 //
260 // Iterators
261 //
262 template <typename PixelT>
264  return _gilView.begin();
265 }
266 
267 template <typename PixelT>
269  return _gilView.end();
270 }
271 
272 template <typename PixelT>
274  return _gilView.rbegin();
275 }
276 
277 template <typename PixelT>
279  return _gilView.rend();
280 }
281 
282 template <typename PixelT>
284  return _gilView.at(x, y);
285 }
286 
287 template <typename PixelT>
289  if (!contiguous) {
290  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Only contiguous == true makes sense");
291  }
292  if (!this->isContiguous()) {
293  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Image's pixels are not contiguous");
294  }
295 
296  return row_begin(0);
297 }
298 
299 template <typename PixelT>
301  if (!contiguous) {
302  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Only contiguous == true makes sense");
303  }
304  if (!this->isContiguous()) {
305  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Image's pixels are not contiguous");
306  }
307 
308  return row_end(getHeight() - 1);
309 }
310 
311 template <typename PixelT>
313  fill_pixels(_gilView, rhs);
314 
315  return *this;
316 }
317 
318 //
319 // On to Image itself. ctors, cctors, and operator=
320 //
321 template <typename PixelT>
322 Image<PixelT>::Image(unsigned int width, unsigned int height, PixelT initialValue)
323  : ImageBase<PixelT>(lsst::geom::ExtentI(width, height)) {
324  *this = initialValue;
325 }
326 
327 template <typename PixelT>
330  *this = initialValue;
331 }
332 
333 template <typename PixelT>
335  *this = initialValue;
336 }
337 
338 template <typename PixelT>
339 Image<PixelT>::Image(Image const& rhs, bool const deep) : ImageBase<PixelT>(rhs, deep) {}
340 // Delegate to copy-constructor for backwards compatibility
341 template <typename PixelT>
342 Image<PixelT>::Image(Image&& rhs) : Image(rhs, false) {}
343 
344 template <typename PixelT>
346  bool const deep)
347  : ImageBase<PixelT>(rhs, bbox, origin, deep) {}
348 
349 template <typename PixelT>
352 
353  return *this;
354 }
355 
356 template <typename PixelT>
359 
360  return *this;
361 }
362 // Delegate to copy-assignment for backwards compatibility
363 template <typename PixelT>
365  return *this = rhs;
366 }
367 
368 #ifndef DOXYGEN // doc for this section has been moved to header
369 
370 template <typename PixelT>
372  lsst::geom::Box2I const& bbox, ImageOrigin origin, bool allowUnsafe) {
373  ImageFitsReader reader(fileName, hdu);
374  *this = reader.read<PixelT>(bbox, origin, allowUnsafe);
375  if (metadata) {
376  metadata->combine(reader.readMetadata());
377  }
378 }
379 
380 template <typename PixelT>
381 Image<PixelT>::Image(fits::MemFileManager& manager, int const hdu,
383  ImageOrigin const origin, bool allowUnsafe) {
384  ImageFitsReader reader(manager, hdu);
385  *this = reader.read<PixelT>(bbox, origin, allowUnsafe);
386  if (metadata) {
387  metadata->combine(reader.readMetadata());
388  }
389 }
390 
391 template <typename PixelT>
393  lsst::geom::Box2I const& bbox, ImageOrigin const origin, bool allowUnsafe) {
394  ImageFitsReader reader(&fitsFile);
395  *this = reader.read<PixelT>(bbox, origin, allowUnsafe);
396  if (metadata) {
397  metadata->combine(reader.readMetadata());
398  }
399 }
400 
401 template <typename PixelT>
402 void Image<PixelT>::writeFits(std::string const& fileName,
404  std::string const& mode) const {
405  fits::Fits fitsfile(fileName, mode, fits::Fits::AUTO_CLOSE | fits::Fits::AUTO_CHECK);
406  writeFits(fitsfile, metadata_i);
407 }
408 
409 template <typename PixelT>
410 void Image<PixelT>::writeFits(fits::MemFileManager& manager,
412  std::string const& mode) const {
413  fits::Fits fitsfile(manager, mode, fits::Fits::AUTO_CLOSE | fits::Fits::AUTO_CHECK);
414  writeFits(fitsfile, metadata_i);
415 }
416 
417 template <typename PixelT>
418 void Image<PixelT>::writeFits(fits::Fits& fitsfile,
420  fitsfile.writeImage(*this, fits::ImageWriteOptions(*this), metadata);
421 }
422 
423 template <typename PixelT>
424 void Image<PixelT>::writeFits(std::string const& filename, fits::ImageWriteOptions const& options,
426  std::shared_ptr<Mask<MaskPixel> const> mask) const {
427  fits::Fits fitsfile(filename, mode, fits::Fits::AUTO_CLOSE | fits::Fits::AUTO_CHECK);
428  writeFits(fitsfile, options, header, mask);
429 }
430 
431 template <typename PixelT>
432 void Image<PixelT>::writeFits(fits::MemFileManager& manager, fits::ImageWriteOptions const& options,
434  std::shared_ptr<Mask<MaskPixel> const> mask) const {
435  fits::Fits fitsfile(manager, mode, fits::Fits::AUTO_CLOSE | fits::Fits::AUTO_CHECK);
436  writeFits(fitsfile, options, header, mask);
437 }
438 
439 template <typename PixelT>
440 void Image<PixelT>::writeFits(fits::Fits& fitsfile, fits::ImageWriteOptions const& options,
442  std::shared_ptr<Mask<MaskPixel> const> mask) const {
443  fitsfile.writeImage(*this, options, header, mask);
444 }
445 
446 #endif // !DOXYGEN
447 
448 template <typename PixelT>
450  using std::swap; // See Meyers, Effective C++, Item 25
452  ; // no private variables to swap
453 }
454 
455 template <typename PixelT>
457  a.swap(b);
458 }
459 
460 // In-place, per-pixel, sqrt().
461 template <typename PixelT>
463  transform_pixels(_getRawView(), _getRawView(),
464  [](PixelT const& l) -> PixelT { return static_cast<PixelT>(std::sqrt(l)); });
465 }
466 
467 template <typename PixelT>
469  transform_pixels(_getRawView(), _getRawView(), [&rhs](PixelT const& l) -> PixelT { return l + rhs; });
470  return *this;
471 }
472 
473 template <typename PixelT>
475  if (this->getDimensions() != rhs.getDimensions()) {
477  (boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
478  this->getHeight() % rhs.getWidth() % rhs.getHeight())
479  .str());
480  }
481  transform_pixels(_getRawView(), rhs._getRawView(), _getRawView(),
482  [](PixelT const& l, PixelT const& r) -> PixelT { return l + r; });
483  return *this;
484 }
485 
486 template <typename PixelT>
488  for (int y = 0; y != this->getHeight(); ++y) {
489  double const yPos = this->indexToPosition(y, Y);
490  double xPos = this->indexToPosition(0, X);
491  for (typename Image<PixelT>::x_iterator ptr = this->row_begin(y), end = this->row_end(y); ptr != end;
492  ++ptr, ++xPos) {
493  *ptr += function(xPos, yPos);
494  }
495  }
496  return *this;
497 }
498 
499 template <typename PixelT>
500 void Image<PixelT>::scaledPlus(double const c, Image<PixelT> const& rhs) {
501  if (this->getDimensions() != rhs.getDimensions()) {
503  (boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
504  this->getHeight() % rhs.getWidth() % rhs.getHeight())
505  .str());
506  }
507  transform_pixels(
508  _getRawView(), rhs._getRawView(), _getRawView(),
509  [&c](PixelT const& l, PixelT const& r) -> PixelT { return l + static_cast<PixelT>(c * r); });
510 }
511 
512 template <typename PixelT>
514  transform_pixels(_getRawView(), _getRawView(), [&rhs](PixelT const& l) -> PixelT { return l - rhs; });
515  return *this;
516 }
517 
518 template <typename PixelT>
520  if (this->getDimensions() != rhs.getDimensions()) {
522  (boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
523  this->getHeight() % rhs.getWidth() % rhs.getHeight())
524  .str());
525  }
526  transform_pixels(_getRawView(), rhs._getRawView(), _getRawView(),
527  [](PixelT const& l, PixelT const& r) -> PixelT { return l - r; });
528  return *this;
529 }
530 
531 template <typename PixelT>
532 void Image<PixelT>::scaledMinus(double const c, Image<PixelT> const& rhs) {
533  if (this->getDimensions() != rhs.getDimensions()) {
535  (boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
536  this->getHeight() % rhs.getWidth() % rhs.getHeight())
537  .str());
538  }
539  transform_pixels(
540  _getRawView(), rhs._getRawView(), _getRawView(),
541  [&c](PixelT const& l, PixelT const& r) -> PixelT { return l - static_cast<PixelT>(c * r); });
542 }
543 
544 template <typename PixelT>
546  for (int y = 0; y != this->getHeight(); ++y) {
547  double const yPos = this->indexToPosition(y, Y);
548  double xPos = this->indexToPosition(0, X);
549  for (typename Image<PixelT>::x_iterator ptr = this->row_begin(y), end = this->row_end(y); ptr != end;
550  ++ptr, ++xPos) {
551  *ptr -= function(xPos, yPos);
552  }
553  }
554  return *this;
555 }
556 
557 template <typename PixelT>
559  transform_pixels(_getRawView(), _getRawView(), [&rhs](PixelT const& l) -> PixelT { return l * rhs; });
560  return *this;
561 }
562 
563 template <typename PixelT>
565  if (this->getDimensions() != rhs.getDimensions()) {
567  (boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
568  this->getHeight() % rhs.getWidth() % rhs.getHeight())
569  .str());
570  }
571  transform_pixels(_getRawView(), rhs._getRawView(), _getRawView(),
572  [](PixelT const& l, PixelT const& r) -> PixelT { return l * r; });
573  return *this;
574 }
575 
576 template <typename PixelT>
577 void Image<PixelT>::scaledMultiplies(double const c, Image<PixelT> const& rhs) {
578  if (this->getDimensions() != rhs.getDimensions()) {
580  (boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
581  this->getHeight() % rhs.getWidth() % rhs.getHeight())
582  .str());
583  }
584  transform_pixels(
585  _getRawView(), rhs._getRawView(), _getRawView(),
586  [&c](PixelT const& l, PixelT const& r) -> PixelT { return l * static_cast<PixelT>(c * r); });
587 }
588 
589 template <typename PixelT>
591  transform_pixels(_getRawView(), _getRawView(), [&rhs](PixelT const& l) -> PixelT { return l / rhs; });
592  return *this;
593 }
594 //
595 // Specialize float and double for efficiency
596 //
597 template <>
599  double const irhs = 1 / rhs;
600  *this *= irhs;
601  return *this;
602 }
603 
604 template <>
606  float const irhs = 1 / rhs;
607  *this *= irhs;
608  return *this;
609 }
610 
611 template <typename PixelT>
613  if (this->getDimensions() != rhs.getDimensions()) {
615  (boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
616  this->getHeight() % rhs.getWidth() % rhs.getHeight())
617  .str());
618  }
619  transform_pixels(_getRawView(), rhs._getRawView(), _getRawView(),
620  [](PixelT const& l, PixelT const& r) -> PixelT { return l / r; });
621  return *this;
622 }
623 
624 template <typename PixelT>
625 void Image<PixelT>::scaledDivides(double const c, Image<PixelT> const& rhs) {
626  if (this->getDimensions() != rhs.getDimensions()) {
628  (boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
629  this->getHeight() % rhs.getWidth() % rhs.getHeight())
630  .str());
631  }
632  transform_pixels(
633  _getRawView(), rhs._getRawView(), _getRawView(),
634  [&c](PixelT const& l, PixelT const& r) -> PixelT { return l / static_cast<PixelT>(c * r); });
635 }
636 
637 namespace {
638 /*
639  * Worker routine for manipulating images;
640  */
641 template <typename LhsPixelT, typename RhsPixelT>
642 struct plusEq : public pixelOp2<LhsPixelT, RhsPixelT> {
643  LhsPixelT operator()(LhsPixelT lhs, RhsPixelT rhs) const override {
644  return static_cast<LhsPixelT>(lhs + rhs);
645  }
646 };
647 
648 template <typename LhsPixelT, typename RhsPixelT>
649 struct minusEq : public pixelOp2<LhsPixelT, RhsPixelT> {
650  LhsPixelT operator()(LhsPixelT lhs, RhsPixelT rhs) const override {
651  return static_cast<LhsPixelT>(lhs - rhs);
652  }
653 };
654 
655 template <typename LhsPixelT, typename RhsPixelT>
656 struct timesEq : public pixelOp2<LhsPixelT, RhsPixelT> {
657  LhsPixelT operator()(LhsPixelT lhs, RhsPixelT rhs) const override {
658  return static_cast<LhsPixelT>(lhs * rhs);
659  }
660 };
661 
662 template <typename LhsPixelT, typename RhsPixelT>
663 struct divideEq : public pixelOp2<LhsPixelT, RhsPixelT> {
664  LhsPixelT operator()(LhsPixelT lhs, RhsPixelT rhs) const override {
665  return static_cast<LhsPixelT>(lhs / rhs);
666  }
667 };
668 } // namespace
669 
670 template <typename LhsPixelT, typename RhsPixelT>
672  for_each_pixel(lhs, rhs, plusEq<LhsPixelT, RhsPixelT>());
673  return lhs;
674 }
675 
676 template <typename LhsPixelT, typename RhsPixelT>
678  for_each_pixel(lhs, rhs, minusEq<LhsPixelT, RhsPixelT>());
679  return lhs;
680 }
681 
682 template <typename LhsPixelT, typename RhsPixelT>
684  for_each_pixel(lhs, rhs, timesEq<LhsPixelT, RhsPixelT>());
685  return lhs;
686 }
687 
688 template <typename LhsPixelT, typename RhsPixelT>
690  for_each_pixel(lhs, rhs, divideEq<LhsPixelT, RhsPixelT>());
691  return lhs;
692 }
693 
696  if (metadata.exists("ZNAXIS1") && metadata.exists("ZNAXIS2")) {
697  dims = lsst::geom::Extent2I(metadata.getAsInt("ZNAXIS1"), metadata.getAsInt("ZNAXIS2"));
698  } else {
699  dims = lsst::geom::Extent2I(metadata.getAsInt("NAXIS1"), metadata.getAsInt("NAXIS2"));
700  }
702  return lsst::geom::Box2I(xy0, dims);
703 }
704 
705 template <typename T1, typename T2>
706 bool imagesOverlap(ImageBase<T1> const& image1, ImageBase<T2> const& image2) {
707 
708  if (image1.getArea() == 0 || image2.getArea() == 0) {
709  // zero-area images cannot overlap.
710  return false;
711  }
712 
713  auto arr1 = image1.getArray();
714  // get the address of the first and one-past-the-last element of arr1 using ndarray iterators;
715  // this works because the iterators for contiguous 1-d ndarray Arrays are just pointers
716  auto beg1Addr = arr1.front().begin();
717  auto end1Addr = arr1.back().end();
718 
719  auto arr2 = image2.getArray();
720  auto beg2Addr = arr2.front().begin();
721  auto end2Addr = arr2.back().end();
722 
723  auto ptrLess = std::less<void const* const>();
724  return ptrLess(beg1Addr, end2Addr) && ptrLess(beg2Addr, end1Addr);
725 }
726 
727 //
728 // Explicit instantiations
729 //
731 #define INSTANTIATE_OPERATOR(OP_EQ, T) \
732  template Image<T>& operator OP_EQ(Image<T>& lhs, Image<std::uint16_t> const& rhs); \
733  template Image<T>& operator OP_EQ(Image<T>& lhs, Image<int> const& rhs); \
734  template Image<T>& operator OP_EQ(Image<T>& lhs, Image<float> const& rhs); \
735  template Image<T>& operator OP_EQ(Image<T>& lhs, Image<double> const& rhs); \
736  template Image<T>& operator OP_EQ(Image<T>& lhs, Image<std::uint64_t> const& rhs);
737 
738 #define INSTANTIATE(T) \
739  template class ImageBase<T>; \
740  template class Image<T>; \
741  INSTANTIATE_OPERATOR(+=, T); \
742  INSTANTIATE_OPERATOR(-=, T); \
743  INSTANTIATE_OPERATOR(*=, T); \
744  INSTANTIATE_OPERATOR(/=, T)
745 
746 #define INSTANTIATE2(T1, T2) template bool imagesOverlap<T1, T2>(ImageBase<T1> const&, ImageBase<T2> const&);
747 
749 INSTANTIATE(int);
750 INSTANTIATE(float);
751 INSTANTIATE(double);
753 
757 INSTANTIATE2(std::uint16_t, double);
759 
761 INSTANTIATE2(int, int);
762 INSTANTIATE2(int, float);
763 INSTANTIATE2(int, double);
765 
767 INSTANTIATE2(float, int);
768 INSTANTIATE2(float, float);
769 INSTANTIATE2(float, double);
771 
772 INSTANTIATE2(double, std::uint16_t);
773 INSTANTIATE2(double, int);
774 INSTANTIATE2(double, float);
775 INSTANTIATE2(double, double);
776 INSTANTIATE2(double, std::uint64_t);
777 
781 INSTANTIATE2(std::uint64_t, double);
783 
785 } // namespace image
786 } // namespace afw
787 } // namespace lsst
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
int end
double x
#define INSTANTIATE(FROMSYS, TOSYS)
Definition: Detector.cc:484
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:49
afw::table::Key< afw::table::Array< MaskPixelT > > mask
#define INSTANTIATE2(ImagePixelT1, ImagePixelT2)
Definition: MaskedImage.cc:695
uint64_t * ptr
Definition: RangeSet.cc:88
int y
Definition: SpanSet.cc:49
table::Key< int > b
table::Key< int > a
A simple struct that combines the two arguments that must be passed to most cfitsio routines and cont...
Definition: fits.h:297
Lifetime-management for memory that goes into FITS memory files.
Definition: fits.h:121
A class used to request that array accesses be checked.
Definition: ImageBase.h:74
std::shared_ptr< daf::base::PropertyList > readMetadata()
Read the image's FITS header.
The base class for all image classed (Image, Mask, MaskedImage, ...)
Definition: ImageBase.h:102
iterator end() const
Return an STL compliant iterator to the end of the image.
Definition: Image.cc:268
iterator begin() const
Return an STL compliant iterator to the start of the image.
Definition: Image.cc:263
static _view_t _allocateView(lsst::geom::Extent2I const &dimensions, Manager::Ptr &manager)
Definition: Image.cc:58
Reference< PixelT >::type PixelReference
A Reference to a PixelT.
Definition: ImageBase.h:117
PixelReference operator()(int x, int y)
Return a reference to the pixel (x, y) in LOCAL coordinates.
Definition: Image.cc:184
static _view_t _makeSubView(lsst::geom::Extent2I const &dimensions, lsst::geom::Extent2I const &offset, const _view_t &view)
Definition: Image.cc:78
int getWidth() const
Return the number of columns in the image.
Definition: ImageBase.h:294
_view_t::reverse_iterator reverse_iterator
An STL compliant reverse iterator.
Definition: ImageBase.h:129
x_iterator fast_iterator
A fast STL compliant iterator for contiguous images N.b.
Definition: ImageBase.h:137
lsst::geom::Box2I getBBox(ImageOrigin origin=PARENT) const
Definition: ImageBase.h:445
ndarray::Array< PixelT, 2, 1 > Array
A mutable ndarray representation of the image.
Definition: ImageBase.h:149
int getArea() const
Return the area of the image.
Definition: ImageBase.h:298
lsst::geom::Extent2I getDimensions() const
Return the image's size; useful for passing to constructors.
Definition: ImageBase.h:356
void assign(ImageBase const &rhs, lsst::geom::Box2I const &bbox=lsst::geom::Box2I(), ImageOrigin origin=PARENT)
Copy pixels from another image to a specified subregion of this image.
Definition: Image.cc:166
int getHeight() const
Return the number of rows in the image.
Definition: ImageBase.h:296
_view_t::iterator iterator
An STL compliant iterator.
Definition: ImageBase.h:125
ImageBase & operator=(const ImageBase &rhs)
Shallow assignment operator.
Definition: Image.cc:153
iterator at(int x, int y) const
Return an STL compliant iterator at the point (x, y)
Definition: Image.cc:283
ConstReference< PixelT >::type PixelConstReference
A ConstReference to a PixelT.
Definition: ImageBase.h:119
reverse_iterator rbegin() const
Return an STL compliant reverse iterator to the start of the image.
Definition: Image.cc:273
_view_t _getRawView() const
Definition: ImageBase.h:465
PixelReference get(lsst::geom::Point2I const &index, ImageOrigin origin)
Return a reference to a single pixel (with no bounds check).
Definition: Image.cc:222
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: ImageBase.h:133
void swap(ImageBase &rhs)
Definition: Image.cc:246
reverse_iterator rend() const
Return an STL compliant reverse iterator to the end of the image.
Definition: Image.cc:278
A FITS reader class for regular Images.
Image< PixelT > read(lsst::geom::Box2I const &bbox=lsst::geom::Box2I(), ImageOrigin origin=PARENT, bool allowUnsafe=false)
Read the Image.
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
void scaledPlus(double const c, Image< PixelT > const &rhs)
Add Image c*rhs to lhs.
Definition: Image.cc:500
Image & operator*=(PixelT const rhs)
Multiply lhs by scalar rhs.
Definition: Image.cc:558
void scaledMinus(double const c, Image< PixelT > const &rhs)
Subtract Image c*rhs from lhs.
Definition: Image.cc:532
Image & operator-=(PixelT const rhs)
Subtract scalar rhs from lhs.
Definition: Image.cc:513
Image & operator=(const PixelT rhs)
Set the image's pixels to rhs.
Definition: Image.cc:350
void scaledMultiplies(double const c, Image< PixelT > const &rhs)
Multiply lhs by Image c*rhs (i.e. pixel-by-pixel multiplication)
Definition: Image.cc:577
Image & operator+=(PixelT const rhs)
Add scalar rhs to lhs.
Definition: Image.cc:468
void swap(Image &rhs)
Definition: Image.cc:449
Image & operator/=(PixelT const rhs)
Divide lhs by scalar rhs.
Definition: Image.cc:590
void writeFits(std::string const &fileName, std::shared_ptr< lsst::daf::base::PropertySet const > metadata=std::shared_ptr< lsst::daf::base::PropertySet const >(), std::string const &mode="w") const
Write an image to a regular FITS file.
void scaledDivides(double const c, Image< PixelT > const &rhs)
Divide lhs by Image c*rhs (i.e. pixel-by-pixel division)
Definition: Image.cc:625
friend class Image
Definition: Image.h:72
A Function taking two arguments.
Definition: Function.h:259
Class for storing generic metadata.
Definition: PropertySet.h:67
int getAsInt(std::string const &name) const
Get the last value for a bool/char/short/int property name (possibly hierarchical).
bool exists(std::string const &name) const
Determine if a name (possibly hierarchical) exists.
An integer coordinate rectangle.
Definition: Box.h:55
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104
Definition: Polygon.cc:25
void swap(CameraSys &a, CameraSys &b)
Definition: CameraSys.h:157
lsst::geom::Point2I getImageXY0FromMetadata(daf::base::PropertySet &metadata, std::string const &wcsName, bool strip=false)
Definition: wcsUtils.cc:98
std::string const wcsNameForXY0
Definition: ImageBase.h:70
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Image< LhsPixelT > & operator+=(Image< LhsPixelT > &lhs, Image< RhsPixelT > const &rhs)
Add lhs to Image rhs (i.e. pixel-by-pixel addition) where types are different.
Definition: Image.cc:671
Image< LhsPixelT > & operator-=(Image< LhsPixelT > &lhs, Image< RhsPixelT > const &rhs)
Subtract lhs from Image rhs (i.e. pixel-by-pixel subtraction) where types are different.
Definition: Image.cc:677
void for_each_pixel(Image< LhsT > &lhs, pixelOp0< LhsT > const &func)
Set each pixel in an Image<LhsT> to func()
lsst::geom::Box2I bboxFromMetadata(daf::base::PropertySet &metadata)
Determine the image bounding box from its metadata (FITS header)
Definition: Image.cc:694
Image< LhsPixelT > & operator/=(Image< LhsPixelT > &lhs, Image< RhsPixelT > const &rhs)
Divide lhs by Image rhs (i.e. pixel-by-pixel division) where types are different.
Definition: Image.cc:689
double indexToPosition(double ind)
Convert image index to image position.
Definition: ImageUtils.h:55
Image< LhsPixelT > & operator*=(Image< LhsPixelT > &lhs, Image< RhsPixelT > const &rhs)
Multiply lhs by Image rhs (i.e. pixel-by-pixel multiplication) where types are different.
Definition: Image.cc:683
bool imagesOverlap(ImageBase< T1 > const &image1, ImageBase< T2 > const &image2)
Return true if the pixels for two images or masks overlap in memory.
Definition: Image.cc:706
void swap(Image< PixelT > &a, Image< PixelT > &b)
Definition: Image.cc:456
Extent< int, 2 > ExtentI
Definition: Extent.h:396
Extent< int, 2 > Extent2I
Definition: Extent.h:397
def writeFits(filename, stamp_ims, metadata, type_name, write_mask, write_variance)
Definition: stamps.py:40
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T sqrt(T... args)
A functor class equivalent to std::function<LhsT (LhsT, RhsT)>, but with a virtual operator()
T swap(T... args)