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
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 
46 #include "lsst/pex/logging/Trace.h"
47 #include "lsst/pex/logging/Log.h"
48 #include "lsst/pex/policy/Policy.h"
50 
51 #include "lsst/ip/diffim.h"
52 
53 namespace afwGeom = lsst::afw::geom;
54 namespace afwImage = lsst::afw::image;
55 namespace afwMath = lsst::afw::math;
56 namespace pexExcept = lsst::pex::exceptions;
57 namespace pexLog = lsst::pex::logging;
58 namespace pexPolicy = lsst::pex::policy;
59 
60 namespace lsst {
61 namespace ip {
62 namespace diffim {
63 
67 template <typename PixelT>
68 Eigen::MatrixXd imageToEigenMatrix(
70  ) {
71  unsigned int rows = img.getHeight();
72  unsigned int cols = img.getWidth();
73  Eigen::MatrixXd M = Eigen::MatrixXd::Zero(rows, cols);
74  for (int y = 0; y != img.getHeight(); ++y) {
75  int x = 0;
76  for (typename afwImage::Image<PixelT>::x_iterator ptr = img.row_begin(y);
77  ptr != img.row_end(y); ++ptr, ++x) {
78  // M is addressed row, col. Need to invert y-axis.
79  // WARNING : CHECK UNIT TESTS BEFORE YOU COMMIT THIS (-y-1) INVERSION
80  M(rows-y-1,x) = *ptr;
81  }
82  }
83  return M;
84 }
85 
86 Eigen::MatrixXi maskToEigenMatrix(
88  ) {
89  unsigned int rows = mask.getHeight();
90  unsigned int cols = mask.getWidth();
91  Eigen::MatrixXi M = Eigen::MatrixXi::Zero(rows, cols);
92  for (int y = 0; y != mask.getHeight(); ++y) {
93  int x = 0;
95  ptr != mask.row_end(y); ++ptr, ++x) {
96  // M is addressed row, col. Need to invert y-axis.
97  // WARNING : CHECK UNIT TESTS BEFORE YOU COMMIT THIS (-y-1) INVERSION
98  M(rows-y-1,x) = *ptr;
99  }
100  }
101  return M;
102 }
103 
104 
120 template <typename PixelT, typename BackgroundT>
122  lsst::afw::image::MaskedImage<PixelT> const &templateImage,
123  lsst::afw::image::MaskedImage<PixelT> const &scienceMaskedImage,
124  lsst::afw::math::Kernel const &convolutionKernel,
125  BackgroundT background,
126  bool invert
127  ) {
128 
129  boost::timer t;
130  t.restart();
131 
132  afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
134  convolutionControl.setDoNormalize(false);
135  afwMath::convolve(convolvedMaskedImage, templateImage,
136  convolutionKernel, convolutionControl);
137 
138  /* Add in background */
139  *(convolvedMaskedImage.getImage()) += background;
140 
141  /* Do actual subtraction */
142  convolvedMaskedImage -= scienceMaskedImage;
143 
144  /* Invert */
145  if (invert) {
146  convolvedMaskedImage *= -1.0;
147  }
148 
149  double time = t.elapsed();
150  pexLog::TTrace<5>("lsst.ip.diffim.convolveAndSubtract",
151  "Total compute time to convolve and subtract : %.2f s", time);
152 
153  return convolvedMaskedImage;
154 }
155 
171 template <typename PixelT, typename BackgroundT>
173  lsst::afw::image::Image<PixelT> const &templateImage,
174  lsst::afw::image::MaskedImage<PixelT> const &scienceMaskedImage,
175  lsst::afw::math::Kernel const &convolutionKernel,
176  BackgroundT background,
177  bool invert
178  ) {
179 
180  boost::timer t;
181  t.restart();
182 
183  afwImage::MaskedImage<PixelT> convolvedMaskedImage(templateImage.getDimensions());
185  convolutionControl.setDoNormalize(false);
186  afwMath::convolve(*convolvedMaskedImage.getImage(), templateImage,
187  convolutionKernel, convolutionControl);
188 
189  /* Add in background */
190  *(convolvedMaskedImage.getImage()) += background;
191 
192  /* Do actual subtraction */
193  *convolvedMaskedImage.getImage() -= *scienceMaskedImage.getImage();
194 
195  /* Invert */
196  if (invert) {
197  *convolvedMaskedImage.getImage() *= -1.0;
198  }
199  *convolvedMaskedImage.getMask() <<= *scienceMaskedImage.getMask();
200  *convolvedMaskedImage.getVariance() <<= *scienceMaskedImage.getVariance();
201 
202  double time = t.elapsed();
203  pexLog::TTrace<5>("lsst.ip.diffim.convolveAndSubtract",
204  "Total compute time to convolve and subtract : %.2f s", time);
205 
206  return convolvedMaskedImage;
207 }
208 
209 /***********************************************************************************************************/
210 //
211 // Explicit instantiations
212 //
213 
214 template
215 Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image<float> const &);
216 
217 template
218 Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image<double> const &);
219 
220 template class FindSetBits<lsst::afw::image::Mask<> >;
221 template class ImageStatistics<float>;
222 template class ImageStatistics<double>;
223 
224 /* */
225 
226 #define p_INSTANTIATE_convolveAndSubtract(TEMPLATE_IMAGE_T, TYPE) \
227  template \
228  lsst::afw::image::MaskedImage<TYPE> convolveAndSubtract( \
229  lsst::afw::image::TEMPLATE_IMAGE_T<TYPE> const& templateImage, \
230  lsst::afw::image::MaskedImage<TYPE> const& scienceMaskedImage, \
231  lsst::afw::math::Kernel const& convolutionKernel, \
232  double background, \
233  bool invert); \
234  \
235  template \
236  afwImage::MaskedImage<TYPE> convolveAndSubtract( \
237  lsst::afw::image::TEMPLATE_IMAGE_T<TYPE> const& templateImage, \
238  lsst::afw::image::MaskedImage<TYPE> const& scienceMaskedImage, \
239  lsst::afw::math::Kernel const& convolutionKernel, \
240  lsst::afw::math::Function2<double> const& backgroundFunction, \
241  bool invert); \
242 
243 #define INSTANTIATE_convolveAndSubtract(TYPE) \
244 p_INSTANTIATE_convolveAndSubtract(Image, TYPE) \
245 p_INSTANTIATE_convolveAndSubtract(MaskedImage, TYPE)
246 /*
247  * Here are the instantiations.
248  *
249  * Do we need double diffim code? It isn't sufficient to remove it here; you'll have to also remove at
250  * least SpatialModelKernel<double> and swig instantiations thereof
251  */
254 
255 }}} // end of namespace lsst::ip::diffim
int y
An include file to include the public header files for lsst::afw::math.
ImagePtr getImage(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:869
An include file to include the header files for lsst::afw::geom.
An include file to include the header files for lsst::ip::diffim.
x_iterator row_begin(int y) const
Definition: Image.h:319
Eigen::MatrixXi maskToEigenMatrix(lsst::afw::image::Mask< lsst::afw::image::MaskPixel > const &mask)
x_iterator row_end(int y) const
Return an x_iterator to the end of the y&#39;th row.
Definition: Image.h:324
Parameters to control convolution.
Definition: ConvolveImage.h:58
definition of the Trace messaging facilities
afwImage::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)
Implement fundamental difference imaging step of convolution and subtraction : D = I - (K*T + bg) whe...
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
An include file to include the header files for lsst::afw::image.
int getWidth() const
Return the number of columns in the image.
Definition: Image.h:237
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:77
int x
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...
MaskPtr getMask(bool const noThrow=false) const
Return a (Ptr to) the MaskedImage&#39;s mask.
Definition: MaskedImage.h:879
#define INSTANTIATE_convolveAndSubtract(TYPE)
int getHeight() const
Return the number of rows in the image.
Definition: Image.h:239
geom::Extent2I getDimensions() const
Definition: MaskedImage.h:904
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:134
geom::Extent2I getDimensions() const
Return the image&#39;s size; useful for passing to constructors.
Definition: Image.h:298
Eigen::MatrixXd imageToEigenMatrix(lsst::afw::image::Image< PixelT > const &img)
Turns Image into a 2-D Eigen Matrix.