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
ImageSubtract.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
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 
34 #include <iostream>
35 #include <numeric>
36 #include <limits>
37 
38 #include "boost/timer.hpp"
39 
40 #include "Eigen/Core"
41 
42 #include "lsst/afw/image.h"
43 #include "lsst/afw/math.h"
44 #include "lsst/log/Log.h"
46 
47 #include "lsst/ip/diffim.h"
48 
49 namespace afwImage = lsst::afw::image;
50 namespace afwMath = lsst::afw::math;
52 namespace dafBase = lsst::daf::base;
53 
54 namespace lsst {
55 namespace ip {
56 namespace diffim {
57 
61 template <typename PixelT>
62 Eigen::MatrixXd imageToEigenMatrix(
64  ) {
65  unsigned int rows = img.getHeight();
66  unsigned int cols = img.getWidth();
67  Eigen::MatrixXd M = Eigen::MatrixXd::Zero(rows, cols);
68  for (int y = 0; y != img.getHeight(); ++y) {
69  int x = 0;
71  ptr != img.row_end(y); ++ptr, ++x) {
72  // M is addressed row, col. Need to invert y-axis.
73  // WARNING : CHECK UNIT TESTS BEFORE YOU COMMIT THIS (-y-1) INVERSION
74  M(rows-y-1,x) = *ptr;
75  }
76  }
77  return M;
78 }
79 
80 Eigen::MatrixXi maskToEigenMatrix(
82  ) {
83  unsigned int rows = mask.getHeight();
84  unsigned int cols = mask.getWidth();
85  Eigen::MatrixXi M = Eigen::MatrixXi::Zero(rows, cols);
86  for (int y = 0; y != mask.getHeight(); ++y) {
87  int x = 0;
89  ptr != mask.row_end(y); ++ptr, ++x) {
90  // M is addressed row, col. Need to invert y-axis.
91  // WARNING : CHECK UNIT TESTS BEFORE YOU COMMIT THIS (-y-1) INVERSION
92  M(rows-y-1,x) = *ptr;
93  }
94  }
95  return M;
96 }
97 
98 
114 template <typename PixelT, typename BackgroundT>
116  lsst::afw::image::MaskedImage<PixelT> const &templateImage,
117  lsst::afw::image::MaskedImage<PixelT> const &scienceMaskedImage,
118  lsst::afw::math::Kernel const &convolutionKernel,
119  BackgroundT background,
120  bool invert
121  ) {
122 
123  boost::timer t;
124  t.restart();
125 
126  afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
128  convolutionControl.setDoNormalize(false);
129  afwMath::convolve(convolvedMaskedImage, templateImage,
130  convolutionKernel, convolutionControl);
131 
132  /* Add in background */
133  *(convolvedMaskedImage.getImage()) += background;
134 
135  /* Do actual subtraction */
136  convolvedMaskedImage -= scienceMaskedImage;
137 
138  /* Invert */
139  if (invert) {
140  convolvedMaskedImage *= -1.0;
141  }
142 
143  double time = t.elapsed();
144  LOGL_DEBUG("TRACE4.ip.diffim.convolveAndSubtract",
145  "Total compute time to convolve and subtract : %.2f s", time);
146 
147  return convolvedMaskedImage;
148 }
149 
165 template <typename PixelT, typename BackgroundT>
167  lsst::afw::image::Image<PixelT> const &templateImage,
168  lsst::afw::image::MaskedImage<PixelT> const &scienceMaskedImage,
169  lsst::afw::math::Kernel const &convolutionKernel,
170  BackgroundT background,
171  bool invert
172  ) {
173 
174  boost::timer t;
175  t.restart();
176 
177  afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
179  convolutionControl.setDoNormalize(false);
180  afwMath::convolve(*convolvedMaskedImage.getImage(), templateImage,
181  convolutionKernel, convolutionControl);
182 
183  /* Add in background */
184  *(convolvedMaskedImage.getImage()) += background;
185 
186  /* Do actual subtraction */
187  *convolvedMaskedImage.getImage() -= *scienceMaskedImage.getImage();
188 
189  /* Invert */
190  if (invert) {
191  *convolvedMaskedImage.getImage() *= -1.0;
192  }
193  convolvedMaskedImage.getMask()->assign(*scienceMaskedImage.getMask());
194  convolvedMaskedImage.getVariance()->assign(*scienceMaskedImage.getVariance());
195 
196  double time = t.elapsed();
197  LOGL_DEBUG("TRACE4.ip.diffim.convolveAndSubtract",
198  "Total compute time to convolve and subtract : %.2f s", time);
199 
200  return convolvedMaskedImage;
201 }
202 
203 /***********************************************************************************************************/
204 //
205 // Explicit instantiations
206 //
207 
208 template
209 Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image<float> const &);
210 
211 template
213 
214 template class FindSetBits<lsst::afw::image::Mask<> >;
215 template class ImageStatistics<float>;
216 template class ImageStatistics<double>;
217 
218 /* */
219 
220 #define p_INSTANTIATE_convolveAndSubtract(TEMPLATE_IMAGE_T, TYPE) \
221  template \
222  lsst::afw::image::MaskedImage<TYPE> convolveAndSubtract( \
223  lsst::afw::image::TEMPLATE_IMAGE_T<TYPE> const& templateImage, \
224  lsst::afw::image::MaskedImage<TYPE> const& scienceMaskedImage, \
225  lsst::afw::math::Kernel const& convolutionKernel, \
226  double background, \
227  bool invert); \
228  \
229  template \
230  afwImage::MaskedImage<TYPE> convolveAndSubtract( \
231  lsst::afw::image::TEMPLATE_IMAGE_T<TYPE> const& templateImage, \
232  lsst::afw::image::MaskedImage<TYPE> const& scienceMaskedImage, \
233  lsst::afw::math::Kernel const& convolutionKernel, \
234  lsst::afw::math::Function2<double> const& backgroundFunction, \
235  bool invert); \
236 
237 #define INSTANTIATE_convolveAndSubtract(TYPE) \
238 p_INSTANTIATE_convolveAndSubtract(Image, TYPE) \
239 p_INSTANTIATE_convolveAndSubtract(MaskedImage, TYPE)
240 /*
241  * Here are the instantiations.
242  *
243  * Do we need double diffim code? It isn't sufficient to remove it here; you'll have to also remove at
244  * least SpatialModelKernel<double>
245  */
248 
249 }}} // end of namespace lsst::ip::diffim
double x
afw::table::Key< afw::table::Array< ImagePixelT > > image
afw::table::Key< afw::table::Array< MaskPixelT > > mask
LSST DM logging module built on log4cxx.
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:504
uint64_t * ptr
Definition: RangeSet.cc:88
int y
Definition: SpanSet.cc:49
int getWidth() const
Return the number of columns in the image.
Definition: ImageBase.h:294
lsst::geom::Extent2I getDimensions() const
Return the image's size; useful for passing to constructors.
Definition: ImageBase.h:356
int getHeight() const
Return the number of rows in the image.
Definition: ImageBase.h:296
x_iterator row_begin(int y) const
Return an x_iterator to the start of the y'th row.
Definition: ImageBase.h:401
x_iterator row_end(int y) const
Return an x_iterator to the end of the y'th row.
Definition: ImageBase.h:404
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
lsst::geom::Extent2I getDimensions() const
Definition: MaskedImage.h:1085
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1079
MaskPtr getMask() const
Return a (shared_ptr to) the MaskedImage's mask.
Definition: MaskedImage.h:1058
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1046
Parameters to control convolution.
Definition: ConvolveImage.h:50
void setDoNormalize(bool doNormalize)
Definition: ConvolveImage.h:66
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:111
Class to accumulate Mask bits.
Definition: FindSetBits.h:53
Class to calculate difference image statistics.
An include file to include the header files for lsst::ip::diffim.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, ConvolutionControl const &convolutionControl=ConvolutionControl())
Convolve an Image or MaskedImage with a Kernel, setting pixels of an existing output image.
Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image< PixelT > const &img)
Turns a 2-d Image into a 2-d Eigen Matrix.
Eigen::MatrixXi maskToEigenMatrix(lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &mask)
lsst::afw::image::MaskedImage< PixelT > convolveAndSubtract(lsst::afw::image::MaskedImage< PixelT > const &templateImage, lsst::afw::image::MaskedImage< PixelT > const &scienceMaskedImage, lsst::afw::math::Kernel const &convolutionKernel, BackgroundT background, bool invert=true)
Execute fundamental task of convolving template and subtracting it from science image.
Relationship invert(Relationship r)
Given the relationship between two sets A and B (i.e.
Definition: Relationship.h:55
A base class for image defects.
#define INSTANTIATE_convolveAndSubtract(TYPE)