// -*- lsst-c++ -*-
/*
* LSST Data Management System
* Copyright 2008-2016 AURA/LSST.
*
* This product includes software developed by the
* LSST Project (http://www.lsst.org/).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the LSST License Statement and
* the GNU General Public License along with this program. If not,
* see .
*/
/**
* \file
* \brief Implementation for ImageBase and Image
*/
#include
#include
#include "boost/mpl/vector.hpp"
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wunused-variable"
#pragma clang diagnostic pop
#include "boost/filesystem/path.hpp"
#include "boost/format.hpp"
#include "boost/gil/gil_all.hpp"
#include "lsst/afw/fits.h"
#include "lsst/afw/image/Image.h"
#include "lsst/afw/image/ImageAlgorithm.h"
#include "lsst/afw/image/Wcs.h"
#include "lsst/afw/image/fits/fits_io.h"
#include "lsst/afw/image/fits/fits_io_mpl.h"
#include "lsst/pex/exceptions.h"
namespace image = lsst::afw::image;
namespace geom = lsst::afw::geom;
/************************************************************************************************************/
template
typename image::ImageBase::_view_t image::ImageBase::_allocateView(
geom::Extent2I const& dimensions, Manager::Ptr& manager) {
if (dimensions.getX() < 0 || dimensions.getY() < 0) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
str(boost::format("Both width and height must be non-negative: %d, %d") %
dimensions.getX() % dimensions.getY()));
}
if (dimensions.getX() != 0 && dimensions.getY() > std::numeric_limits::max() / dimensions.getX()) {
throw LSST_EXCEPT(pex::exceptions::LengthError,
str(boost::format("Image dimensions (%d x %d) too large; "
"int overflow detected.") %
dimensions.getX() % dimensions.getY()));
}
std::pair r =
ndarray::SimpleManager::allocate(dimensions.getX() * dimensions.getY());
manager = r.first;
return boost::gil::interleaved_view(dimensions.getX(), dimensions.getY(),
(typename _view_t::value_type*)r.second,
dimensions.getX() * sizeof(PixelT));
}
template
typename image::ImageBase::_view_t image::ImageBase::_makeSubView(
geom::Extent2I const& dimensions, geom::Extent2I const& offset, const _view_t& view) {
if (offset.getX() < 0 || offset.getY() < 0 || offset.getX() + dimensions.getX() > view.width() ||
offset.getY() + dimensions.getY() > view.height()) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Box2I(Point2I(%d,%d),Extent2I(%d,%d)) "
"doesn't fit in image %dx%d") %
offset.getX() % offset.getY() % dimensions.getX() % dimensions.getY() %
view.width() % view.height())
.str());
}
return boost::gil::subimage_view(view, offset.getX(), offset.getY(), dimensions.getX(),
dimensions.getY());
}
/**
* Allocator Constructor
*
* allocate a new image with the specified dimensions.
* Sets origin at (0,0)
*/
template
image::ImageBase::ImageBase(geom::Extent2I const& dimensions)
: lsst::daf::base::Citizen(typeid(this)),
_origin(0, 0),
_manager(),
_gilView(_allocateView(dimensions, _manager)) {}
/**
* Allocator Constructor
*
* allocate a new image with the specified dimensions and origin
*/
template
image::ImageBase::ImageBase(geom::Box2I const& bbox)
: lsst::daf::base::Citizen(typeid(this)),
_origin(bbox.getMin()),
_manager(),
_gilView(_allocateView(bbox.getDimensions(), _manager)) {}
/**
* Copy constructor.
*
* \note Unless \c deep is \c true, the new %image will share the old %image's
* pixels;
* this may not be what you want. See also assign(rhs) to copy pixels between
* Image%s
*/
template
image::ImageBase::ImageBase(ImageBase const& rhs, ///< Right-hand-side %image
bool const deep ///< If false, new ImageBase shares storage with rhs;
///< if true make a new, standalone, ImageBase
)
: lsst::daf::base::Citizen(typeid(this)),
_origin(rhs._origin),
_manager(rhs._manager),
_gilView(rhs._gilView) {
if (deep) {
ImageBase tmp(getBBox());
tmp.assign(*this); // now copy the pixels
swap(tmp);
}
}
/**
* Copy constructor to make a copy of part of an %image.
*
* The bbox ignores X0/Y0 if origin == LOCAL, and uses it if origin == PARENT.
*
* \note Unless \c deep is \c true, the new %image will share the old %image's
* pixels;
* this is probably what you want
*/
template
image::ImageBase::ImageBase(ImageBase const& rhs, ///< Right-hand-side %image
geom::Box2I const& bbox, ///< Specify desired region
ImageOrigin const origin, ///< Specify the coordinate system of the bbox
bool const deep ///< If false, new ImageBase shares storage with rhs;
///< if true make a new, standalone, ImageBase
)
: lsst::daf::base::Citizen(typeid(this)),
_origin((origin == PARENT) ? bbox.getMin() : rhs._origin + geom::Extent2I(bbox.getMin())),
_manager(rhs._manager), // reference counted pointer, don't copy pixels
_gilView(_makeSubView(bbox.getDimensions(), _origin - rhs._origin, rhs._gilView)) {
if (deep) {
ImageBase tmp(getBBox());
tmp.assign(*this); // now copy the pixels
swap(tmp);
}
}
/**
* Construction from ndarray::Array and NumPy.
*
* \note ndarray and NumPy indexes are ordered (y,x), but Image indices are
* ordered (x,y).
*
* Unless deep is true, the new image will share memory with the array if the
* the
* dimension is contiguous in memory. If the last dimension is not contiguous,
* the array
* will be deep-copied in Python, but the constructor will fail to compile in
* pure C++.
*/
template
image::ImageBase::ImageBase(Array const& array, bool deep, geom::Point2I const& xy0)
: lsst::daf::base::Citizen(typeid(this)),
_origin(xy0),
_manager(array.getManager()),
_gilView(boost::gil::interleaved_view(array.template getSize<1>(), array.template getSize<0>(),
(typename _view_t::value_type*)array.getData(),
array.template getStride<0>() * sizeof(PixelT))) {
if (deep) {
ImageBase tmp(*this, true);
swap(tmp);
}
}
/// Shallow assignment operator.
///
/// \note that this has the effect of making the lhs share pixels with the rhs
/// which may
/// not be what you intended; to copy the pixels, use assign(rhs)
///
/// \note this behaviour is required to make the swig interface work, otherwise
/// I'd
/// declare this function private
template
image::ImageBase& image::ImageBase::operator=(ImageBase const& rhs) {
ImageBase tmp(rhs);
swap(tmp); // See Meyers, Effective C++, Item 11
return *this;
}
/// Set the lhs's %pixel values to equal the rhs's
///
/// \deprecated use assign(rhs) instead
template
image::ImageBase& image::ImageBase::operator<<=(ImageBase const& rhs) {
assign(rhs);
return *this;
}
/**
* Copy pixels from another image to a specified subregion of this image.
*
* \param[in] rhs source image whose pixels are to be copied into this image
* (the destination)
* \param[in] bbox subregion of this image to set; if empty (the default) then
* all pixels are set
* \param[in] origin origin of bbox: if PARENT then the lower left pixel of
* this image is at xy0
* if LOCAL then the lower left pixel of this image is at 0,0
*
* \throw lsst::pex::exceptions::LengthError if the dimensions of rhs and the
* specified subregion of
* this image do not match.
*/
template
void image::ImageBase::assign(ImageBase const& rhs, geom::Box2I const& bbox, ImageOrigin origin) {
auto lhsDim = bbox.isEmpty() ? getDimensions() : bbox.getDimensions();
if (lhsDim != rhs.getDimensions()) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Dimension mismatch: %dx%d v. %dx%d") % lhsDim.getX() %
lhsDim.getY() % rhs.getWidth() % rhs.getHeight())
.str());
}
if (bbox.isEmpty()) {
copy_pixels(rhs._gilView, _gilView);
} else {
auto lhsOff = (origin == PARENT) ? bbox.getMin() - _origin : geom::Extent2I(bbox.getMin());
auto lhsGilView = _makeSubView(lhsDim, lhsOff, _gilView);
copy_pixels(rhs._gilView, lhsGilView);
}
}
/// Return a reference to the pixel (x, y)
template
typename image::ImageBase::PixelReference image::ImageBase::operator()(int x, int y) {
return const_cast::PixelReference>(
static_cast::PixelConstReference>(_gilView(x, y)[0]));
}
/// Return a reference to the pixel (x, y) with bounds checking
template
typename image::ImageBase::PixelReference image::ImageBase::operator()(
int x, int y, image::CheckIndices const& check) {
if (check && (x < 0 || x >= getWidth() || y < 0 || y >= getHeight())) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Index (%d, %d) is out of range [0--%d], [0--%d]") % x % y %
(getWidth() - 1) % (getHeight() - 1))
.str());
}
return const_cast::PixelReference>(
static_cast::PixelConstReference>(_gilView(x, y)[0]));
}
/// Return a const reference to the pixel (x, y)
template
typename image::ImageBase::PixelConstReference image::ImageBase::operator()(int x,
int y) const {
return _gilView(x, y)[0];
}
/// Return a const reference to the pixel (x, y) with bounds checking
template
typename image::ImageBase::PixelConstReference image::ImageBase::operator()(
int x, int y, image::CheckIndices const& check) const {
if (check && (x < 0 || x >= getWidth() || y < 0 || y >= getHeight())) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Index (%d, %d) is out of range [0--%d], [0--%d]") % x % y %
(this->getWidth() - 1) % (this->getHeight() - 1))
.str());
}
return _gilView(x, y)[0];
}
template
void image::ImageBase::swap(ImageBase& rhs) {
using std::swap; // See Meyers, Effective C++, Item 25
swap(_manager, rhs._manager); // just swapping the pointers
swap(_gilView, rhs._gilView);
swap(_origin, rhs._origin);
}
template
void image::swap(ImageBase& a, ImageBase& b) {
a.swap(b);
}
template
typename image::ImageBase::Array image::ImageBase::getArray() {
int rowStride = reinterpret_cast(row_begin(1)) - reinterpret_cast(row_begin(0));
return ndarray::external(reinterpret_cast(row_begin(0)),
ndarray::makeVector(getHeight(), getWidth()), ndarray::makeVector(rowStride, 1),
this->_manager);
}
template
typename image::ImageBase::ConstArray image::ImageBase::getArray() const {
int rowStride = reinterpret_cast(row_begin(1)) - reinterpret_cast(row_begin(0));
return ndarray::external(reinterpret_cast(row_begin(0)),
ndarray::makeVector(getHeight(), getWidth()), ndarray::makeVector(rowStride, 1),
this->_manager);
}
//
// Iterators
//
/// Return an STL compliant iterator to the start of the %image
///
/// Note that this isn't especially efficient; see \link imageIterators\endlink
/// for
/// a discussion
template
typename image::ImageBase::iterator image::ImageBase::begin() const {
return _gilView.begin();
}
/// Return an STL compliant iterator to the end of the %image
template
typename image::ImageBase::iterator image::ImageBase::end() const {
return _gilView.end();
}
/// Return an STL compliant reverse iterator to the start of the %image
template
typename image::ImageBase::reverse_iterator image::ImageBase::rbegin() const {
return _gilView.rbegin();
}
/// Return an STL compliant reverse iterator to the end of the %image
template
typename image::ImageBase::reverse_iterator image::ImageBase::rend() const {
return _gilView.rend();
}
/// Return an STL compliant iterator at the point (x, y)
template
typename image::ImageBase::iterator image::ImageBase::at(int x, int y) const {
return _gilView.at(x, y);
}
/// Return a fast STL compliant iterator to the start of the %image which must
/// be contiguous
///
/// \exception lsst::pex::exceptions::Runtime
/// Argument \a contiguous is false, or the pixels are not in fact contiguous
template
typename image::ImageBase::fast_iterator image::ImageBase::begin(
bool contiguous ///< Pixels are contiguous (must be true)
) const {
if (!contiguous) {
throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Only contiguous == true makes sense");
}
if (!this->isContiguous()) {
throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Image's pixels are not contiguous");
}
return row_begin(0);
}
/// Return a fast STL compliant iterator to the end of the %image which must be
/// contiguous
///
/// \exception lsst::pex::exceptions::Runtime
/// Argument \a contiguous is false, or the pixels are not in fact contiguous
template
typename image::ImageBase::fast_iterator image::ImageBase::end(
bool contiguous ///< Pixels are contiguous (must be true)
) const {
if (!contiguous) {
throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Only contiguous == true makes sense");
}
if (!this->isContiguous()) {
throw LSST_EXCEPT(lsst::pex::exceptions::RuntimeError, "Image's pixels are not contiguous");
}
return row_end(getHeight() - 1);
}
/************************************************************************************************************/
/// Set the %image's pixels to rhs
template
image::ImageBase& image::ImageBase::operator=(PixelT const rhs) {
fill_pixels(_gilView, rhs);
return *this;
}
/************************************************************************************************************/
//
// On to Image itself. ctors, cctors, and operator=
//
/**
* Create an initialised Image of the specified size
*
* \note Many lsst::afw::image and lsst::afw::math objects define a \c
* dimensions member
* which may be conveniently used to make objects of an appropriate size
*/
template
image::Image::Image(unsigned int width, ///< number of columns
unsigned int height, ///< number of rows
PixelT initialValue ///< Initial value
)
: image::ImageBase(geom::ExtentI(width, height)) {
*this = initialValue;
}
/**
* Create an initialised Image of the specified size
*
* \note Many lsst::afw::image and lsst::afw::math objects define a \c
* dimensions member
* which may be conveniently used to make objects of an appropriate size
*/
template
image::Image::Image(geom::Extent2I const& dimensions, ///< Number of columns, rows
PixelT initialValue ///< Initial value
)
: image::ImageBase(dimensions) {
*this = initialValue;
}
/**
* Create an initialized Image of the specified size
*/
template
image::Image::Image(geom::Box2I const& bbox, ///< dimensions and origin of desired Image
PixelT initialValue ///< Initial value
)
: image::ImageBase(bbox) {
*this = initialValue;
}
/**
* Copy constructor.
*
* \note Unless \c deep is \c true, the new %image will share the old %image's
* pixels;
* this may not be what you want. See also assign(rhs) to copy pixels between
* Image%s
*/
template
image::Image::Image(Image const& rhs, ///< Right-hand-side Image
bool const deep ///< If false, new Image shares storage with rhs; if true
///< make a new, standalone, ImageBase
)
: image::ImageBase(rhs, deep) {}
/**
* Copy constructor to make a copy of part of an Image.
*
* The bbox ignores X0/Y0 if origin == LOCAL, and uses it if origin == PARENT.
*
* \note Unless \c deep is \c true, the new %image will share the old %image's
* pixels;
* this is probably what you want
*/
template
image::Image::Image(Image const& rhs, ///< Right-hand-side Image
geom::Box2I const& bbox, ///< Specify desired region
ImageOrigin const origin, ///< Coordinate system of the bbox
bool const deep ///< If false, new ImageBase shares storage with rhs; if true
///< make a new, standalone, ImageBase
)
: image::ImageBase(rhs, bbox, origin, deep) {}
/// Set the %image's pixels to rhs
template
image::Image& image::Image::operator=(PixelT const rhs) {
this->ImageBase::operator=(rhs);
return *this;
}
/// Assignment operator.
///
/// \note that this has the effect of making the lhs share pixels with the rhs
/// which may
/// not be what you intended; to copy the pixels, use assign(rhs)
///
/// \note this behaviour is required to make the swig interface work, otherwise
/// I'd
/// declare this function private
template
image::Image& image::Image::operator=(Image const& rhs) {
this->ImageBase::operator=(rhs);
return *this;
}
/************************************************************************************************************/
#ifndef DOXYGEN // doc for this section has been moved to header
template
image::Image::Image(std::string const& fileName, int hdu, PTR(daf::base::PropertySet) metadata,
geom::Box2I const& bbox, ImageOrigin origin)
: image::ImageBase() {
fits::Fits fitsfile(fileName, "r", fits::Fits::AUTO_CLOSE | fits::Fits::AUTO_CHECK);
fitsfile.setHdu(hdu);
try {
*this = Image(fitsfile, metadata, bbox, origin);
} catch (lsst::afw::fits::FitsError& e) {
fitsfile.status = 0; // reset so we can read NAXIS
if (fitsfile.getImageDim() == 0) { // no pixels to read
LSST_EXCEPT_ADD(e, str(boost::format("HDU %d has NAXIS == 0") % hdu));
}
throw e;
}
}
template
image::Image::Image(fits::MemFileManager& manager, int const hdu,
PTR(daf::base::PropertySet) metadata, geom::Box2I const& bbox,
ImageOrigin const origin)
: image::ImageBase() {
fits::Fits fitsfile(manager, "r", fits::Fits::AUTO_CLOSE | fits::Fits::AUTO_CHECK);
fitsfile.setHdu(hdu);
*this = Image(fitsfile, metadata, bbox, origin);
}
template
image::Image::Image(fits::Fits& fitsfile, PTR(daf::base::PropertySet) metadata,
geom::Box2I const& bbox, ImageOrigin const origin)
: image::ImageBase() {
typedef boost::mpl::vector
fits_image_types;
if (!metadata) {
metadata.reset(new daf::base::PropertyList());
}
fits_read_image(fitsfile, *this, *metadata, bbox, origin);
}
template
void image::Image::writeFits(std::string const& fileName,
CONST_PTR(lsst::daf::base::PropertySet) metadata_i,
std::string const& mode) const {
fits::Fits fitsfile(fileName, mode, fits::Fits::AUTO_CLOSE | fits::Fits::AUTO_CHECK);
writeFits(fitsfile, metadata_i);
}
template
void image::Image::writeFits(fits::MemFileManager& manager,
CONST_PTR(lsst::daf::base::PropertySet) metadata_i,
std::string const& mode) const {
fits::Fits fitsfile(manager, mode, fits::Fits::AUTO_CLOSE | fits::Fits::AUTO_CHECK);
writeFits(fitsfile, metadata_i);
}
template
void image::Image::writeFits(fits::Fits& fitsfile,
CONST_PTR(lsst::daf::base::PropertySet) metadata_i) const {
PTR(daf::base::PropertySet) metadata;
PTR(daf::base::PropertySet)
wcsAMetadata = image::detail::createTrivialWcsAsPropertySet(image::detail::wcsNameForXY0, this->getX0(),
this->getY0());
if (metadata_i) {
metadata = metadata_i->deepCopy();
metadata->combine(wcsAMetadata);
} else {
metadata = wcsAMetadata;
}
image::fits_write_image(fitsfile, *this, metadata);
}
#endif // !DOXYGEN
/************************************************************************************************************/
template
void image::Image::swap(Image& rhs) {
using std::swap; // See Meyers, Effective C++, Item 25
ImageBase::swap(rhs);
; // no private variables to swap
}
template
void image::swap(Image& a, Image& b) {
a.swap(b);
}
/************************************************************************************************************/
// In-place, per-pixel, sqrt().
template
void image::Image::sqrt() {
transform_pixels(_getRawView(), _getRawView(),
[](PixelT const& l) -> PixelT { return static_cast(std::sqrt(l)); });
}
/// Add scalar rhs to lhs
template
image::Image& image::Image::operator+=(PixelT const rhs) {
transform_pixels(_getRawView(), _getRawView(), [&rhs](PixelT const& l) -> PixelT { return l + rhs; });
return *this;
}
/// Add Image rhs to lhs
template
image::Image& image::Image::operator+=(Image const& rhs) {
if (this->getDimensions() != rhs.getDimensions()) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
this->getHeight() % rhs.getWidth() % rhs.getHeight())
.str());
}
transform_pixels(_getRawView(), rhs._getRawView(), _getRawView(),
[](PixelT const& l, PixelT const& r) -> PixelT { return l + r; });
return *this;
}
/**
* @brief Add a Function2(x, y) to an Image
*/
template
image::Image& image::Image::operator+=(
lsst::afw::math::Function2 const& function ///< function to add
) {
for (int y = 0; y != this->getHeight(); ++y) {
double const yPos = this->indexToPosition(y, image::Y);
double xPos = this->indexToPosition(0, image::X);
for (typename Image::x_iterator ptr = this->row_begin(y), end = this->row_end(y); ptr != end;
++ptr, ++xPos) {
*ptr += function(xPos, yPos);
}
}
return *this;
}
/// Add Image c*rhs to lhs
template
void image::Image::scaledPlus(double const c, Image const& rhs) {
if (this->getDimensions() != rhs.getDimensions()) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
this->getHeight() % rhs.getWidth() % rhs.getHeight())
.str());
}
transform_pixels(
_getRawView(), rhs._getRawView(), _getRawView(),
[&c](PixelT const& l, PixelT const& r) -> PixelT { return l + static_cast(c * r); });
}
/// Subtract scalar rhs from lhs
template
image::Image& image::Image::operator-=(PixelT const rhs) {
transform_pixels(_getRawView(), _getRawView(), [&rhs](PixelT const& l) -> PixelT { return l - rhs; });
return *this;
}
/// Subtract Image rhs from lhs
template
image::Image& image::Image::operator-=(Image const& rhs) {
if (this->getDimensions() != rhs.getDimensions()) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
this->getHeight() % rhs.getWidth() % rhs.getHeight())
.str());
}
transform_pixels(_getRawView(), rhs._getRawView(), _getRawView(),
[](PixelT const& l, PixelT const& r) -> PixelT { return l - r; });
return *this;
}
/// Subtract Image c*rhs from lhs
template
void image::Image::scaledMinus(double const c, Image const& rhs) {
if (this->getDimensions() != rhs.getDimensions()) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
this->getHeight() % rhs.getWidth() % rhs.getHeight())
.str());
}
transform_pixels(
_getRawView(), rhs._getRawView(), _getRawView(),
[&c](PixelT const& l, PixelT const& r) -> PixelT { return l - static_cast(c * r); });
}
/**
* @brief Subtract a Function2(x, y) from an Image
*/
template
image::Image& image::Image::operator-=(
lsst::afw::math::Function2 const& function ///< function to add
) {
for (int y = 0; y != this->getHeight(); ++y) {
double const yPos = this->indexToPosition(y, image::Y);
double xPos = this->indexToPosition(0, image::X);
for (typename Image::x_iterator ptr = this->row_begin(y), end = this->row_end(y); ptr != end;
++ptr, ++xPos) {
*ptr -= function(xPos, yPos);
}
}
return *this;
}
/// Multiply lhs by scalar rhs
template
image::Image& image::Image::operator*=(PixelT const rhs) {
transform_pixels(_getRawView(), _getRawView(), [&rhs](PixelT const& l) -> PixelT { return l * rhs; });
return *this;
}
/// Multiply lhs by Image rhs (i.e. %pixel-by-%pixel multiplication)
template
image::Image& image::Image::operator*=(Image const& rhs) {
if (this->getDimensions() != rhs.getDimensions()) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
this->getHeight() % rhs.getWidth() % rhs.getHeight())
.str());
}
transform_pixels(_getRawView(), rhs._getRawView(), _getRawView(),
[](PixelT const& l, PixelT const& r) -> PixelT { return l * r; });
return *this;
}
/// Multiply lhs by Image c*rhs (i.e. %pixel-by-%pixel multiplication)
template
void image::Image::scaledMultiplies(double const c, Image const& rhs) {
if (this->getDimensions() != rhs.getDimensions()) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
this->getHeight() % rhs.getWidth() % rhs.getHeight())
.str());
}
transform_pixels(
_getRawView(), rhs._getRawView(), _getRawView(),
[&c](PixelT const& l, PixelT const& r) -> PixelT { return l * static_cast(c * r); });
}
/// Divide lhs by scalar rhs
///
/// \note Floating point types implement this by multiplying by the 1/rhs
template
image::Image& image::Image::operator/=(PixelT const rhs) {
transform_pixels(_getRawView(), _getRawView(), [&rhs](PixelT const& l) -> PixelT { return l / rhs; });
return *this;
}
//
// Specialize float and double for efficiency
//
namespace lsst {
namespace afw {
namespace image {
template <>
Image& Image::operator/=(double const rhs) {
double const irhs = 1 / rhs;
*this *= irhs;
return *this;
}
template <>
Image& Image::operator/=(float const rhs) {
float const irhs = 1 / rhs;
*this *= irhs;
return *this;
}
}
}
}
/// Divide lhs by Image rhs (i.e. %pixel-by-%pixel division)
template
image::Image& image::Image::operator/=(Image const& rhs) {
if (this->getDimensions() != rhs.getDimensions()) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
this->getHeight() % rhs.getWidth() % rhs.getHeight())
.str());
}
transform_pixels(_getRawView(), rhs._getRawView(), _getRawView(),
[](PixelT const& l, PixelT const& r) -> PixelT { return l / r; });
return *this;
}
/// Divide lhs by Image c*rhs (i.e. %pixel-by-%pixel division)
template
void image::Image::scaledDivides(double const c, Image const& rhs) {
if (this->getDimensions() != rhs.getDimensions()) {
throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
(boost::format("Images are of different size, %dx%d v %dx%d") % this->getWidth() %
this->getHeight() % rhs.getWidth() % rhs.getHeight())
.str());
}
transform_pixels(
_getRawView(), rhs._getRawView(), _getRawView(),
[&c](PixelT const& l, PixelT const& r) -> PixelT { return l / static_cast(c * r); });
}
/************************************************************************************************************/
namespace {
/*
* Worker routine for manipulating images;
*/
template
struct plusEq : public lsst::afw::image::pixelOp2 {
LhsPixelT operator()(LhsPixelT lhs, RhsPixelT rhs) const { return static_cast(lhs + rhs); }
};
template
struct minusEq : public lsst::afw::image::pixelOp2 {
LhsPixelT operator()(LhsPixelT lhs, RhsPixelT rhs) const { return static_cast(lhs - rhs); }
};
template
struct timesEq : public lsst::afw::image::pixelOp2 {
LhsPixelT operator()(LhsPixelT lhs, RhsPixelT rhs) const { return static_cast(lhs * rhs); }
};
template
struct divideEq : public lsst::afw::image::pixelOp2 {
LhsPixelT operator()(LhsPixelT lhs, RhsPixelT rhs) const { return static_cast(lhs / rhs); }
};
}
/// Add lhs to Image rhs (i.e. %pixel-by-%pixel addition) where types are
/// different
///
template
image::Image& image::operator+=(image::Image& lhs, image::Image const& rhs) {
image::for_each_pixel(lhs, rhs, plusEq());
return lhs;
}
/// Subtract lhs from Image rhs (i.e. %pixel-by-%pixel subtraction) where types
/// are different
///
template
image::Image& image::operator-=(image::Image& lhs, image::Image const& rhs) {
image::for_each_pixel(lhs, rhs, minusEq());
return lhs;
}
/// Multiply lhs by Image rhs (i.e. %pixel-by-%pixel multiplication) where types
/// are different
///
template
image::Image& image::operator*=(image::Image& lhs, image::Image const& rhs) {
image::for_each_pixel(lhs, rhs, timesEq());
return lhs;
}
/// Divide lhs by Image rhs (i.e. %pixel-by-%pixel division) where types are
/// different
///
template
image::Image& image::operator/=(image::Image& lhs, image::Image const& rhs) {
image::for_each_pixel(lhs, rhs, divideEq());
return lhs;
}
/************************************************************************************************************/
//
// Explicit instantiations
//
/// \cond
#define INSTANTIATE_OPERATOR(OP_EQ, T) \
template image::Image& image::operator OP_EQ(image::Image& lhs, \
image::Image const& rhs); \
template image::Image& image::operator OP_EQ(image::Image& lhs, image::Image const& rhs); \
template image::Image& image::operator OP_EQ(image::Image& lhs, image::Image const& rhs); \
template image::Image& image::operator OP_EQ(image::Image& lhs, image::Image const& rhs); \
template image::Image& image::operator OP_EQ(image::Image& lhs, \
image::Image const& rhs);
#define INSTANTIATE(T) \
template class image::ImageBase; \
template class image::Image; \
INSTANTIATE_OPERATOR(+=, T); \
INSTANTIATE_OPERATOR(-=, T); \
INSTANTIATE_OPERATOR(*=, T); \
INSTANTIATE_OPERATOR(/=, T)
INSTANTIATE(std::uint16_t);
INSTANTIATE(int);
INSTANTIATE(float);
INSTANTIATE(double);
INSTANTIATE(std::uint64_t);
/// \endcond