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
Shapelet.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 
25 #include "Eigen/Core"
26 
27 #include "lsst/afw/geom.h"
28 #include "lsst/afw/image/Wcs.h"
37 
38 namespace lsst {
39 namespace meas {
40 namespace algorithms {
41 
42  // All of the functionality is imported from shapelet::BVec
43  // Just repeat the constructors, destructors, and op=
45  {
46  public :
47 
51 
52  ShapeletImpl(int order, double sigma) :
53  BVec(order,sigma), _cov()
54  {}
55 
56  ShapeletImpl(int order, double sigma, const ShapeletVector& vector) :
57  BVec(order,sigma,vector), _cov()
58  {}
59 
60  ShapeletImpl(int order, double sigma, const ShapeletVector& vector, const ShapeletCovariance& cov) :
61  BVec(order,sigma,vector), _cov(new ShapeletCovariance(cov))
62  {}
63 
64  ShapeletImpl(const BVec& rhs) :
65  BVec(rhs), _cov()
66  {}
67 
68  ShapeletImpl(const ShapeletImpl& rhs) :
69  BVec(rhs)
70  {
71  if (rhs._cov.get()) _cov.reset(new ShapeletCovariance(*rhs._cov));
72  }
73 
75 
76  using BVec::size;
77 
78  void operator=(const ShapeletImpl& rhs)
79  {
80  BVec::operator=(rhs);
81  if (rhs._cov.get()) {
82  if (_cov.get() &&
83  _cov->TMV_colsize() == rhs._cov->TMV_colsize() &&
84  _cov->TMV_rowsize() == rhs._cov->TMV_rowsize()) {
85  *_cov = *rhs._cov;
86  } else {
87  _cov.reset(new ShapeletCovariance(*rhs._cov));
88  }
89  } else {
90  _cov.reset();
91  }
92  }
93 
94  boost::shared_ptr<ShapeletCovariance>& getCovariance()
95  { return _cov; }
96 
97  std::complex<double> getPQ(int p, int q)
98  {
99  using shapelet::DVector;
100 
101  if (p < q) return std::conj(getPQ(q,p));
102  else {
103  const DVector& v = BVec::getValues();
104  if (p == q) {
105  int k = p*(2*p+3);
106  return v(k);
107  } else {
108  int k = (p+q)*(p+q+1)/2 + 2*q;
109  return std::complex<double>(v(k),v(k+1));
110  }
111  }
112  }
113 
114  double evaluateAt(double x, double y)
115  {
116  // TODO: This is not efficient, but the functionality
117  // is required for Kernel.
118  //
119  // It is much more efficient to run makePsi on many
120  // positions at once. Furthermore, makePsi is optimized
121  // to be run on many positions at the same time, which
122  // involved some coding choices that probably make it even
123  // less efficient for a single position.
124  //
125  // However, if this is a significant time component,
126  // then a more efficient version that is intended for
127  // a single position could be written.
128  using shapelet::CDVector;
129  using shapelet::DMatrix;
130  using shapelet::makePsi;
131 
132  std::complex<double> z(x,y);
133  z /= getSigma();
134  DMatrix psi(1,BVec::size());
135  CDVector zList(1);
136  zList(0) = z;
137  makePsi(psi,TMV_vview(zList),BVec::getOrder(),0);
138  return (psi * getValues())(0);
139  }
140 
141  private :
142  boost::shared_ptr<ShapeletCovariance> _cov;
143  };
144 
145  Shapelet::Shapelet(int order, double sigma) :
146  pImpl(new ShapeletImpl(order,sigma))
147  {}
148 
149  Shapelet::Shapelet(int order, double sigma, const ShapeletVector& vector) :
150  pImpl(new ShapeletImpl(order,sigma,vector))
151  {}
152 
153  Shapelet::Shapelet(int order, double sigma, const ShapeletVector& vector, const ShapeletCovariance& cov) :
154  pImpl(new ShapeletImpl(order,sigma,vector,cov))
155  {}
156 
158  pImpl(new ShapeletImpl(*rhs.pImpl))
159  {}
160 
162  { delete pImpl; pImpl = 0; }
163 
165  {
166  *pImpl = *rhs.pImpl;
167  return *this;
168  }
169 
170  int Shapelet::getOrder() const
171  { return pImpl->getOrder(); }
172 
173  double Shapelet::getSigma() const
174  { return pImpl->getSigma(); }
175 
176  int Shapelet::size() const
177  { return pImpl->size(); }
178 
180  { return pImpl->getValues(); }
181 
183  { return pImpl->getCovariance().get(); }
184 
185  boost::shared_ptr<const Shapelet::ShapeletCovariance> Shapelet::getCovariance() const
186  { return pImpl->getCovariance(); }
187 
189  { return pImpl->setSigma(sigma); }
190 
191  std::complex<double> Shapelet::getPQ(int p, int q)
192  { return pImpl->getPQ(p,q); }
193 
194  double Shapelet::evaluateAt(const PointD& pos)
195  { return evaluateAt(pos.getX(),pos.getY()); }
196 
197  double Shapelet::evaluateAt(double x, double y)
198  {
199  // TODO: MJ - I couldn't reproduce the crash with g++ 4.3.3.
200  // Need to try this on lsst computers I guess.
201 #if 0 // crashes g++ 4.3.5
202  return pImpl->evaluateAt(x,y);
203 #else
204  // MJ - This is definitely wrong! It starts an infinite loop!
205  return evaluateAt(lsst::afw::geom::PointD(x, y));
206 #endif
207  }
208 
209  Eigen::Matrix2d getJacobian(
210  const lsst::afw::image::Wcs& wcs,
211  const lsst::afw::geom::PointD& pos)
212  {
213  lsst::afw::geom::AffineTransform localTransform =
215 
216  // J = ( du/dx du/dy )
217  // ( dv/dx dv/dy )
218  // where (u,v) are sky coordinates, and (x,y) are chip coordinates.
219  Eigen::Matrix2d J;
220  // Answer comes out in degrees for u,v. *3600 to get arcsec.
221  J(0,0) = localTransform.getMatrix()(0,0) * 3600.;
222  J(0,1) = localTransform.getMatrix()(0,1) * 3600.;
223  J(1,0) = localTransform.getMatrix()(1,0) * 3600.;
224  J(1,1) = localTransform.getMatrix()(1,1) * 3600.;
225  return J;
226  }
227 
228  // Can't use the version in Pixel.h, since we need to use the
229  // LSST Image and Wcs objects, rather than my Image and Trasnsformation.
230  static void getPixList(
231  shapelet::PixelList& pix,
232  const Shapelet::Source& source,
233  const Shapelet::PointD& cen,
234  double aperture,
235  const Shapelet::Exposure& exposure,
237  {
238  using shapelet::Pixel;
239  using shapelet::PixelList;
241 
242  Shapelet::Exposure::MaskedImageT const maskedImage = exposure.getMaskedImage();
243  Shapelet::Exposure::MaskedImageT::Image::ConstPtr imagePtr = maskedImage.getImage();
244  Shapelet::Exposure::MaskedImageT::Mask::ConstPtr maskPtr = maskedImage.getMask();
245  Shapelet::Exposure::MaskedImageT::Variance::ConstPtr variancePtr = maskedImage.getVariance();
246 
247  PointD pos(source.getX(),source.getY());
248  Eigen::Matrix2d J = getJacobian(*(exposure.getWcs()), pos);
249 
250  double det = std::abs(J.determinant());
251  double pixScale = sqrt(det); // arcsec / pixel
252  xdbg<<"pixscale = "<<pixScale<<std::endl;
253 
254  // xAp,yAp are the maximum deviation from the center in x,y
255  // such that u^2+v^2 = aperture^2
256  double xAp = aperture / det *
257  sqrt(J(0,0)*J(0,0) + J(0,1)*J(0,1));
258  double yAp = aperture / det *
259  sqrt(J(1,0)*J(1,0) + J(1,1)*J(1,1));
260  xdbg<<"aperture = "<<aperture<<std::endl;
261  xdbg<<"xap = "<<xAp<<", yap = "<<yAp<<std::endl;
262 
263  double xCen = cen.getX();
264  double yCen = cen.getY();
265  xdbg<<"cen = "<<xCen<<" "<<yCen<<std::endl;
266 
267  // Find the square range bounding the aperture
268  int i1 = int(floor(xCen-xAp));
269  int i2 = int(ceil(xCen+xAp));
270  int j1 = int(floor(yCen-yAp));
271  int j2 = int(ceil(yCen+yAp));
272  xdbg<<"i1,i2,j1,j2 = "<<i1<<','<<i2<<','<<j1<<','<<j2<<std::endl;
273 
274  // Keep the range used withing the borders of the image.
275  // TODO: Need to check what the LSST conventions are for x,y in
276  // image(x,y).
277  // Are they always relative to whatever the lower left position is?
278  // Or relative to the same (0,0) that makes getX0 and getY0
279  // possibly non-zero.
280  int xMin = exposure.getX0();
281  int yMin = exposure.getY0();
282  int xMax = xMin + exposure.getWidth(); // no image->getX1() method?
283  int yMax = yMin + exposure.getHeight();
284  xdbg<<"xMin, yMin = "<<xMin<<" "<<yMin<<std::endl;
285  if (i1 < xMin) { i1 = xMin; }
286  if (i2 > xMax) { i2 = xMax; }
287  if (j1 < yMin) { j1 = yMin; }
288  if (j2 > yMax) { j2 = yMax; }
289  xdbg<<"i1,i2,j1,j2 => "<<i1<<','<<i2<<','<<j1<<','<<j2<<std::endl;
290 
291  // Do this next loop in two passes. First figure out which
292  // pixels we want to use. Then we can resize pix to the full size
293  // we will need, and go back through and enter the pixels.
294  // This saves us a lot of resizing calls in vector, which are
295  // both slow and can fragment the memory.
296  xdbg<<"nx = "<<i2-i1+1<<std::endl;
297  xdbg<<"ny = "<<j2-j1+1<<std::endl;
298  Assert(i2-i1+1 >= 0);
299  Assert(j2-j1+1 >= 0);
300  std::vector<std::vector<bool> > shouldUsePix(
301  i2-i1+1,std::vector<bool>(j2-j1+1,false));
302  int nPix = 0;
303 
304  const double apsq = aperture*aperture;
305 
306  double chipX = i1-xCen;
307  for(int i=i1;i<=i2;++i,chipX+=1.) {
308  double chipY = j1-yCen;
309  double u = J(0,0)*chipX+J(0,1)*chipY;
310  double v = J(1,0)*chipX+J(1,1)*chipY;
311  for(int j=j1;j<=j2;++j,u+=J(0,1),v+=J(1,1)) {
312  // u,v are in arcsec
313  double rsq = u*u + v*v;
314  if ( ((*maskPtr)(i,j) & ~okmask) &&
315  (rsq <= apsq) ) {
316  shouldUsePix[i-i1][j-j1] = true;
317  ++nPix;
318  }
319  }
320  }
321 
322  xdbg<<"npix = "<<nPix<<std::endl;
323  pix.resize(nPix);
324 
325  xdbg<<"pixlist size = "<<nPix<<" = "<<nPix*sizeof(Pixel)<<" bytes = "
326  <<nPix*sizeof(Pixel)/1024.<<" KB\n";
327 
328  // Now the real loop that stores the flux values.
329  double sky = 0.0; // FIXME: this should be source.getSky() (maybe?)
330 
331  int k=0;
332  chipX = i1-xCen;
333  for(int i=i1;i<=i2;++i,chipX+=1.) {
334  double chipY = j1-yCen;
335  double u = J(0,0)*chipX+J(0,1)*chipY;
336  double v = J(1,0)*chipX+J(1,1)*chipY;
337  for(int j=j1;j<=j2;++j,u+=J(0,1),v+=J(1,1)) {
338  if (shouldUsePix[i-i1][j-j1]) {
339  double flux = (*imagePtr)(i,j)-sky;
340  double variance = (*variancePtr)(i,j);
341  if (variance > 0.0) {
342  double inverseSigma = sqrt(1.0/variance);
343  Assert(k < int(pix.size()));
344  pix[k++] = Pixel(u,v,flux,inverseSigma);
345  }
346  }
347  }
348  }
349  Assert(k <= int(pix.size()));
350  // Not necessarily k == pix.size()
351  // because we skip pixels with 0.0 variance
352  pix.resize(k); // clear off the extras, if any.
353  Assert(k == int(pix.size()));
354  nPix = pix.size(); // may have changed.
355  xdbg<<"npix => "<<nPix<<std::endl;
356  }
357 
359  const Source& source, const PointD& pos,
360  bool isCentroidFixed, bool isSigmaFixed, double aperture,
361  const Exposure& exposure,
362  const lsst::afw::image::MaskPixel okmask)
363  {
364  using shapelet::Ellipse;
365  using shapelet::PixelList;
366 
367  std::vector<PixelList> pix(1);
368  // Fill PixelList with pixel data around position pos:
369  getPixList(pix[0], source, pos, aperture, exposure, okmask);
370 
371  double sigma = pImpl->getSigma();
372  Ellipse ell;
373  ell.fixGam();
374  if (isCentroidFixed) ell.fixCen();
375  else {
376  // Initial crude estimates to get close to the right value
377  // in case we have a poor starting point.
378  // TODO: These might not be necessary for LSST.
379  // Should compare speed with and without this step.
380  ell.peakCentroid(pix[0],aperture/3.);
381  ell.crudeMeasure(pix[0],sigma);
382  }
383  if (isSigmaFixed) ell.fixMu();
384  long flag = 0;
385  int order = getOrder();
386  if (!isCentroidFixed || !isSigmaFixed) {
387  if (!ell.measure(pix,order,order+4,order,sigma,flag,1.e-4)) {
388  return false;
389  }
390  if (flag) return false;
391  if (!isSigmaFixed) {
392  double mu = ell.getMu();
393  sigma *= exp(mu);
394  dbg<<"sigma = "<<sigma<<std::endl;
395  Assert(sigma > 0);
396  pImpl->setSigma(sigma);
397  }
398  }
399 
400  ShapeletCovariance* cov = pImpl->getCovariance().get();
401  if (ell.measureShapelet(pix,*pImpl,order,order+4,order,cov)) {
402  return true;
403  } else {
404  return false;
405  }
406  }
407 
408  Shapelet::Shapelet(const shapelet::BVec& bvec) : pImpl(new ShapeletImpl(bvec)) {}
409 
410  const shapelet::BVec& Shapelet::viewAsBVec() const { return *pImpl; }
412 
413 }}}
414 
int y
Matrix const getMatrix() const
boost::uint16_t MaskPixel
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.
const ShapeletVector & getValues() const
get the values as a vector.
Definition: Shapelet.cc:179
~Shapelet()
Destructor needs to delete pImpl.
Definition: Shapelet.cc:161
void setSigma(double sigma)
set a new value of sigma
Definition: Shapelet.cc:188
geom::AffineTransform linearizePixelToSky(coord::Coord const &coord, geom::AngleUnit skyUnit=geom::degrees) const
Return the local linear approximation to Wcs::pixelToSky at a point given in sky coordinates.
Definition: Wcs.cc:934
#define TMV_vview(v)
Definition: MyMatrix.h:202
BVec & operator=(const AssignableToBVec &rhs)
Definition: BVec.cc:39
A class to contain the data, WCS, and other information needed to describe an image of the sky...
Definition: Exposure.h:48
double evaluateAt(const PointD &pos)
Evaluate f(x,y)
Definition: Shapelet.cc:194
Defines the Shapelet class.
Shapelet::ShapeletCovariance ShapeletCovariance
Definition: Shapelet.cc:49
void operator=(const ShapeletImpl &rhs)
Definition: Shapelet.cc:78
tbl::Key< int > wcs
bool measureFromImage(const Source &source, const PointD &pos, bool isCentroidFixed, bool isSigmaFixed, double aperture, const Exposure &exposure, const MaskPixel okmask=0)
measure shapelet decomposition of an image
Definition: Shapelet.cc:358
lsst::afw::geom::PointD PointD
Definition: Shapelet.h:102
Implementation of the WCS standard for a any projection.
Definition: Wcs.h:107
double getSigma() const
get the scale size of the shapelet
Definition: Shapelet.cc:173
int getOrder() const
get the order of the shapelet
Definition: Shapelet.cc:170
const shapelet::BVec & viewAsBVec() const
View the shapelet as a BVec.
Definition: Shapelet.cc:410
Extent< int, N > ceil(Extent< double, N > const &input)
Shapelet & operator=(const Shapelet &rhs)
op= does a deep copy
Definition: Shapelet.cc:164
ShapeletImpl(int order, double sigma)
Definition: Shapelet.cc:52
Eigen::MatrixXd ShapeletCovariance
Definition: Shapelet.h:98
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
boost::shared_ptr< ShapeletCovariance > _cov
Definition: Shapelet.cc:142
void makePsi(DMatrix &psi, CDVectorView z, int order, const DVector *coeff=0)
Definition: PsiHelper.cc:226
std::complex< double > getPQ(int p, int q)
Definition: Shapelet.cc:97
lsst::afw::table::SourceRecord Source
Definition: Shapelet.h:100
ShapeletImpl(const ShapeletImpl &rhs)
Definition: Shapelet.cc:68
AngleUnit const degrees
Definition: Angle.h:92
bool hasCovariance() const
does the shapelet have a covariance matrix stored?
Definition: Shapelet.cc:182
boost::shared_ptr< const ShapeletCovariance > getCovariance() const
get the covariance matrix
Definition: Shapelet.cc:185
Eigen::Matrix2d getJacobian(const lsst::afw::image::Wcs &wcs, const lsst::afw::geom::PointD &pos)
a helper function to deal with the fact that Wcs doesn&#39;t directly return a Jacobian.
Definition: Shapelet.cc:209
An affine coordinate transformation consisting of a linear transformation and an offset.
void setSigma(double sigma)
Definition: BVec.h:69
std::complex< double > getPQ(int p, int q)
Get a complex b_pq value.
Definition: Shapelet.cc:191
lsst::afw::image::Exposure< PixelT > Exposure
Definition: Shapelet.h:101
double evaluateAt(double x, double y)
Definition: Shapelet.cc:114
const DVector & getValues() const
Definition: BVec.h:66
int x
Shapelet::ShapeletVector ShapeletVector
Definition: Shapelet.cc:48
Extent< int, N > floor(Extent< double, N > const &input)
MaskedImage< ImageT, MaskT, VarianceT > MaskedImageT
Definition: Exposure.h:51
ShapeletImpl(int order, double sigma, const ShapeletVector &vector)
Definition: Shapelet.cc:56
Point< double, 2 > PointD
Definition: Point.h:285
Eigen::VectorXd ShapeletVector
Definition: Shapelet.h:97
Record class that contains measurements made on a single exposure.
Definition: Source.h:81
Implementation of the Class MaskedImage.
#define xdbg
Definition: dbg.h:70
ShapeletImpl(int order, double sigma, const ShapeletVector &vector, const ShapeletCovariance &cov)
Definition: Shapelet.cc:60
#define dbg
Definition: dbg.h:69
int size() const
the size of the shapelet vector
Definition: Shapelet.cc:176
boost::shared_ptr< ShapeletCovariance > & getCovariance()
Definition: Shapelet.cc:94
#define Assert(x)
Definition: dbg.h:73
Shapelet(int order, double sigma)
Basic constructor requires order and sigma.
Definition: Shapelet.cc:145