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
ShapeletKernel.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 #include "assert.h"
25 
29 #include "lsst/afw/math/Function.h"
30 
31 namespace lsst {
32 namespace meas {
33 namespace algorithms {
34 
36  public lsst::afw::math::Function2<double>
37  {
38  public :
40 
41  typedef boost::shared_ptr<ShapeletKernelFunction> Ptr;
42  typedef boost::shared_ptr<const ShapeletKernelFunction> ConstPtr;
43 
45  base(shapelet->size()+1),
46  _order(shapelet->getOrder()), _size(shapelet->size())
47  {
48  const Shapelet::ShapeletVector& vals = shapelet->getValues();
49  std::vector<double> v(_size+1);
50  for(int i=0;i<_size;++i) v[i] = vals(i);
51  v[_size] = shapelet->getSigma();
53  }
54 
55  ShapeletKernelFunction(int order, int size) :
56  base(size+1), _order(order), _size(size)
57  {}
58 
59  base::Ptr clone() const
60  { return base::Ptr(new ShapeletKernelFunction(*this)); }
61 
62  double operator()(double x, double y) const
63  {
65  const std::vector<double>& v = base::getParameters();
66  for(int i=0;i<_size;++i) vals(i) = v[i];
67  double sigma = v[_size];
68  Shapelet shapelet(_order,sigma,vals);
69  return shapelet.evaluateAt(x,y);
70  }
71 
72  private :
73  const int _order;
74  const int _size;
75  };
76 
78  public lsst::afw::math::Function2<double>
79  {
80  public :
82 
83  typedef boost::shared_ptr<ShapeletSpatialFunction> Ptr;
84  typedef boost::shared_ptr<const ShapeletSpatialFunction> ConstPtr;
85 
88  ) :
89  base(interp->getFitSize()), _interp(interp), _i(i)
90  {}
91 
92  base::Ptr clone() const
93  { return base::Ptr(new ShapeletSpatialFunction(*this)); }
94 
95  double operator()(double x, double y) const
96  { return _interp->interpolateSingleElement(x,y,_i); }
97 
98  private :
100  const int _i;
101 
102  };
103 
104  // Use a radius of 5 sigma for width,height if they are not specified.
105  template <class ShapeletPtr>
106  inline int getImageSize(ShapeletPtr shapelet, int size)
107  { return size == 0 ? int(ceil(shapelet->getSigma()*10.)) : size; }
108 
110  Shapelet::ConstPtr shapelet, const Wcs::ConstPtr& wcsPtr, const Extent& size) :
111  base( getImageSize(shapelet,size.getX()),
112  getImageSize(shapelet,size.getY()),
113  ShapeletKernelFunction(shapelet) ),
114  _shapelet(shapelet), _wcsPtr(wcsPtr->clone())
115  {}
116 
118  Shapelet::ConstPtr shapelet, const Wcs::ConstPtr& wcsPtr) :
119  base( getImageSize(shapelet,0),
120  getImageSize(shapelet,0),
121  ShapeletKernelFunction(shapelet) ),
122  _shapelet(shapelet), _wcsPtr(wcsPtr->clone())
123  {}
124 
126  Image& image, bool doNormalize, double /*x*/, double /*y*/) const
127  {
128  //TODO This isn't quite right. I need to account for the WCS.
129  using shapelet::CDVector;
130  using shapelet::DMatrix;
131  using shapelet::DVector;
132  using shapelet::makePsi;
133 
134  const int xCen = base::getCtrX();
135  const int yCen = base::getCtrY();
136  const int nX = base::getWidth();
137  const int nY = base::getHeight();
138  const int nPixels = nX * nY;
139  const double sigma = _shapelet->getSigma();
140  CDVector zList(nPixels);
141  int k=0;
142  for(int i=0;i<nX;++i) {
143  for(int j=0;j<nY;++j) {
144  std::complex<double> z(i-xCen,j-yCen);
145  z /= sigma;
146  zList[k++] = z;
147  }
148  }
149  Assert(k == nPixels);
150  DMatrix psi(nPixels,_shapelet->size());
151  makePsi(psi,TMV_vview(zList),_shapelet->getOrder());
152  DVector flux = psi * _shapelet->getValues();
153 
154  double sum = flux.TMV_sumElements();
155  if (doNormalize) flux /= sum;
156 
157  k=0;
158  for(int i=0;i<nX;++i) {
159  for(int j=0;j<nY;++j) {
160  image(i,j) = flux(k++);
161  }
162  }
163 
164  return sum;
165  }
166 
167  std::vector<lsst::afw::math::Kernel::SpatialFunctionPtr> buildSetOfSpatialFunctions(
169  {
170  const int nParam = interp->getFitSize();
172  std::vector<SFPtr> ret;
173  ret.reserve(nParam);
174  for(int i=0;i<nParam;++i) {
175  ret.push_back(SFPtr(new ShapeletSpatialFunction(interp,i)));
176  }
177  return ret;
178  }
179 
182  const Wcs::ConstPtr& wcsPtr, const Extent& size
183  ) :
184  base( getImageSize(interp,size.getX()),
185  getImageSize(interp,size.getY()),
186  ShapeletKernelFunction(interp->getOrder(),interp->getSize()),
187  buildSetOfSpatialFunctions(interp) ),
188  _interp(interp), _wcsPtr(wcsPtr->clone())
189  {}
190 
192  ShapeletInterpolation::ConstPtr interp, const Wcs::ConstPtr& wcsPtr
193  ) :
194  base( getImageSize(interp,0),
195  getImageSize(interp,0),
196  ShapeletKernelFunction(interp->getOrder(),interp->getSize()),
197  buildSetOfSpatialFunctions(interp) ),
198  _interp(interp), _wcsPtr(wcsPtr->clone())
199  {}
200 
202  const Point& pos) const
203  {
204  Shapelet::ConstPtr localShapelet(_interp->interpolate(pos));
205  LocalShapeletKernel::ConstPtr localShapeletKernel(
206  new LocalShapeletKernel(localShapelet, _wcsPtr));
207  return localShapeletKernel;
208  }
209 
211  Image& image, bool doNormalize, double x, double y) const
212  {
213  Point point(x,y);
214  return getLocalKernel(point)->computeImage(image,doNormalize);
215  }
216 
217 }}}
int y
double computeImage(Image &image, bool doNormalize, double x=0.0, double y=0.0) const
Make an image of the kernel at a specified location.
std::vector< lsst::afw::math::Kernel::SpatialFunctionPtr > buildSetOfSpatialFunctions(ShapeletInterpolation::ConstPtr interp)
double computeImage(Image &image, bool doNormalize, double x=0.0, double y=0.0) const
Make an image of the kernel.
LocalShapeletKernel(Shapelet::ConstPtr shapelet, const Wcs::ConstPtr &wcsPtr, const Extent &size)
Constructor from a Shapelet.
int getWidth() const
Return the Kernel&#39;s width.
Definition: Kernel.h:240
boost::shared_ptr< ShapeletSpatialFunction > Ptr
base::Ptr clone() const
Return a pointer to a deep copy of this function.
boost::shared_ptr< lsst::afw::math::Function2< double > > SpatialFunctionPtr
Definition: Kernel.h:143
std::vector< double > const & getParameters() const
Return all function parameters.
Definition: Function.h:145
#define TMV_vview(v)
Definition: MyMatrix.h:202
lsst::afw::math::Function2< double > base
ShapeletInterpolation::ConstPtr _interp
int getHeight() const
Return the Kernel&#39;s height.
Definition: Kernel.h:247
boost::shared_ptr< const LocalShapeletKernel > ConstPtr
ShapeletKernelFunction(Shapelet::ConstPtr shapelet)
double evaluateAt(const PointD &pos)
Evaluate f(x,y)
Definition: Shapelet.cc:194
const ShapeletInterpolation::ConstPtr _interp
Extent< int, N > ceil(Extent< double, N > const &input)
A Function taking two arguments.
Definition: Function.h:300
boost::shared_ptr< const Shapelet > ConstPtr
Definition: Shapelet.h:95
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
Definition: operators.h:1250
int getImageSize(ShapeletPtr shapelet, int size)
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
void makePsi(DMatrix &psi, CDVectorView z, int order, const DVector *coeff=0)
Definition: PsiHelper.cc:226
boost::shared_ptr< Wcs const > ConstPtr
Definition: Wcs.h:114
boost::shared_ptr< const ShapeletKernelFunction > ConstPtr
double operator()(double x, double y) const
boost::shared_ptr< ShapeletKernelFunction > Ptr
boost::shared_ptr< const ShapeletSpatialFunction > ConstPtr
ShapeletSpatialFunction(ShapeletInterpolation::ConstPtr interp, int i)
int x
base::Ptr clone() const
Return a pointer to a deep copy of this function.
boost::shared_ptr< const ShapeletInterpolation > ConstPtr
LocalShapeletKernel::ConstPtr getLocalKernel(const Point &pos) const
Get the LocalShapeletKernel at a given point.
lsst::afw::math::Function2< double > base
Define the basic Function classes.
Eigen::VectorXd ShapeletVector
Definition: Shapelet.h:97
ShapeletKernel(ShapeletInterpolation::ConstPtr interp, const Wcs::ConstPtr &wcsPtr, const Extent &size)
Constructor from a ShapeletInterpolation.
boost::shared_ptr< Function2< double > > Ptr
Definition: Function.h:304
double operator()(double x, double y) const
A kernel described by a function.
Definition: Kernel.h:628
int getCtrY() const
Return y index of kernel&#39;s center.
Definition: Kernel.h:272
A class to represent a 2-dimensional array of pixels.
Definition: PSF.h:43
Defines LocalShapeletKernel and ShapeletKernel.
#define Assert(x)
Definition: dbg.h:73
void setParameters(std::vector< double > const &params)
Set all function parameters.
Definition: Function.h:175
int getCtrX() const
Return x index of kernel&#39;s center.
Definition: Kernel.h:263