LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
LSSTDataManagementBasePackage
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/afw/geom.h"
45 #include "lsst/log/Log.h"
46 #include "lsst/pex/policy/Policy.h"
48 
49 #include "lsst/ip/diffim.h"
50 
51 namespace afwGeom = lsst::afw::geom;
52 namespace afwImage = lsst::afw::image;
53 namespace afwMath = lsst::afw::math;
54 namespace pexExcept = lsst::pex::exceptions;
55 namespace pexPolicy = lsst::pex::policy;
56 
57 namespace lsst {
58 namespace ip {
59 namespace diffim {
60 
64 template <typename PixelT>
65 Eigen::MatrixXd imageToEigenMatrix(
67  ) {
68  unsigned int rows = img.getHeight();
69  unsigned int cols = img.getWidth();
70  Eigen::MatrixXd M = Eigen::MatrixXd::Zero(rows, cols);
71  for (int y = 0; y != img.getHeight(); ++y) {
72  int x = 0;
73  for (typename afwImage::Image<PixelT>::x_iterator ptr = img.row_begin(y);
74  ptr != img.row_end(y); ++ptr, ++x) {
75  // M is addressed row, col. Need to invert y-axis.
76  // WARNING : CHECK UNIT TESTS BEFORE YOU COMMIT THIS (-y-1) INVERSION
77  M(rows-y-1,x) = *ptr;
78  }
79  }
80  return M;
81 }
82 
83 Eigen::MatrixXi maskToEigenMatrix(
85  ) {
86  unsigned int rows = mask.getHeight();
87  unsigned int cols = mask.getWidth();
88  Eigen::MatrixXi M = Eigen::MatrixXi::Zero(rows, cols);
89  for (int y = 0; y != mask.getHeight(); ++y) {
90  int x = 0;
92  ptr != mask.row_end(y); ++ptr, ++x) {
93  // M is addressed row, col. Need to invert y-axis.
94  // WARNING : CHECK UNIT TESTS BEFORE YOU COMMIT THIS (-y-1) INVERSION
95  M(rows-y-1,x) = *ptr;
96  }
97  }
98  return M;
99 }
100 
101 
117 template <typename PixelT, typename BackgroundT>
119  lsst::afw::image::MaskedImage<PixelT> const &templateImage,
120  lsst::afw::image::MaskedImage<PixelT> const &scienceMaskedImage,
121  lsst::afw::math::Kernel const &convolutionKernel,
122  BackgroundT background,
123  bool invert
124  ) {
125 
126  boost::timer t;
127  t.restart();
128 
129  afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
131  convolutionControl.setDoNormalize(false);
132  afwMath::convolve(convolvedMaskedImage, templateImage,
133  convolutionKernel, convolutionControl);
134 
135  /* Add in background */
136  *(convolvedMaskedImage.getImage()) += background;
137 
138  /* Do actual subtraction */
139  convolvedMaskedImage -= scienceMaskedImage;
140 
141  /* Invert */
142  if (invert) {
143  convolvedMaskedImage *= -1.0;
144  }
145 
146  double time = t.elapsed();
147  LOGL_DEBUG("TRACE4.ip.diffim.convolveAndSubtract",
148  "Total compute time to convolve and subtract : %.2f s", time);
149 
150  return convolvedMaskedImage;
151 }
152 
168 template <typename PixelT, typename BackgroundT>
170  lsst::afw::image::Image<PixelT> const &templateImage,
171  lsst::afw::image::MaskedImage<PixelT> const &scienceMaskedImage,
172  lsst::afw::math::Kernel const &convolutionKernel,
173  BackgroundT background,
174  bool invert
175  ) {
176 
177  boost::timer t;
178  t.restart();
179 
180  afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
182  convolutionControl.setDoNormalize(false);
183  afwMath::convolve(*convolvedMaskedImage.getImage(), templateImage,
184  convolutionKernel, convolutionControl);
185 
186  /* Add in background */
187  *(convolvedMaskedImage.getImage()) += background;
188 
189  /* Do actual subtraction */
190  *convolvedMaskedImage.getImage() -= *scienceMaskedImage.getImage();
191 
192  /* Invert */
193  if (invert) {
194  *convolvedMaskedImage.getImage() *= -1.0;
195  }
196  convolvedMaskedImage.getMask()->assign(*scienceMaskedImage.getMask());
197  convolvedMaskedImage.getVariance()->assign(*scienceMaskedImage.getVariance());
198 
199  double time = t.elapsed();
200  LOGL_DEBUG("TRACE4.ip.diffim.convolveAndSubtract",
201  "Total compute time to convolve and subtract : %.2f s", time);
202 
203  return convolvedMaskedImage;
204 }
205 
206 /***********************************************************************************************************/
207 //
208 // Explicit instantiations
209 //
210 
211 template
212 Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image<float> const &);
213 
214 template
215 Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image<double> const &);
216 
217 template class FindSetBits<lsst::afw::image::Mask<> >;
218 template class ImageStatistics<float>;
219 template class ImageStatistics<double>;
220 
221 /* */
222 
223 #define p_INSTANTIATE_convolveAndSubtract(TEMPLATE_IMAGE_T, TYPE) \
224  template \
225  lsst::afw::image::MaskedImage<TYPE> convolveAndSubtract( \
226  lsst::afw::image::TEMPLATE_IMAGE_T<TYPE> const& templateImage, \
227  lsst::afw::image::MaskedImage<TYPE> const& scienceMaskedImage, \
228  lsst::afw::math::Kernel const& convolutionKernel, \
229  double background, \
230  bool invert); \
231  \
232  template \
233  afwImage::MaskedImage<TYPE> convolveAndSubtract( \
234  lsst::afw::image::TEMPLATE_IMAGE_T<TYPE> const& templateImage, \
235  lsst::afw::image::MaskedImage<TYPE> const& scienceMaskedImage, \
236  lsst::afw::math::Kernel const& convolutionKernel, \
237  lsst::afw::math::Function2<double> const& backgroundFunction, \
238  bool invert); \
239 
240 #define INSTANTIATE_convolveAndSubtract(TYPE) \
241 p_INSTANTIATE_convolveAndSubtract(Image, TYPE) \
242 p_INSTANTIATE_convolveAndSubtract(MaskedImage, TYPE)
243 /*
244  * Here are the instantiations.
245  *
246  * Do we need double diffim code? It isn't sufficient to remove it here; you'll have to also remove at
247  * least SpatialModelKernel<double>
248  */
251 
252 }}} // end of namespace lsst::ip::diffim
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage&#39;s variance.
Definition: MaskedImage.h:1091
uint64_t * ptr
Definition: RangeSet.cc:88
An include file to include the header files for lsst::ip::diffim.
Eigen::MatrixXi maskToEigenMatrix(lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &mask)
int getHeight() const
Return the number of rows in the image.
Definition: ImageBase.h:316
Parameters to control convolution.
Definition: ConvolveImage.h:50
int y
Definition: SpanSet.cc:49
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
x_iterator row_begin(int y) const
Return an x_iterator to the start of the y&#39;th row.
Definition: ImageBase.h:419
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:489
Class to accumulate Mask bits.
Definition: FindSetBits.h:53
lsst::geom::Extent2I getDimensions() const
Definition: MaskedImage.h:1097
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:1058
LSST DM logging module built on log4cxx.
A base class for image defects.
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:78
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:74
Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image< PixelT > const &img)
Turns a 2-d Image into a 2-d Eigen Matrix.
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...
#define INSTANTIATE_convolveAndSubtract(TYPE)
double x
Class to calculate difference image statistics.
MaskPtr getMask() const
Return a (shared_ptr to) the MaskedImage&#39;s mask.
Definition: MaskedImage.h:1070
afw::table::Key< afw::table::Array< MaskPixelT > > mask
int getWidth() const
Return the number of columns in the image.
Definition: ImageBase.h:314
afw::table::Key< afw::table::Array< ImagePixelT > > image
x_iterator row_end(int y) const
Return an x_iterator to the end of the y&#39;th row.
Definition: ImageBase.h:422
lsst::geom::Extent2I getDimensions() const
Return the image&#39;s size; useful for passing to constructors.
Definition: ImageBase.h:374
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:112
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:59