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
CrudeMeasure.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 <vector>
26 
30 
31 namespace lsst {
32 namespace meas {
33 namespace algorithms {
34 namespace shapelet {
35 
36  class CrudeSolver : public NLSolver
37  {
38  public :
40  const PixelList& pix, double sigma, double I1,
41  DVector& xinit);
43 
44  void calculateF(const DVector& x, DVector& f) const;
45  void calculateJ(
46  const DVector& x, const DVector& f,
47  DMatrix& df) const;
48 
49  private :
50  double _sigma;
54  mutable CDVector _Z1;
55  mutable DVector _E;
56  mutable DVector _Rsq;
57  mutable DVector _f1;
58  double _I1;
60  };
61 
63  const PixelList& pix,
64  double sigma, double I1, DVector& xInit) :
65  _sigma(sigma), _I(pix.size()), _Z(pix.size()), _W(pix.size()),
66  _Z1(pix.size()), _E(pix.size()), _Rsq(pix.size()), _f1(pix.size()),
67  _I1(I1), _xInit(xInit)
68  {
69  const int nPix = pix.size();
70  for(int i=0;i<nPix;++i) {
71  _Z(i) = pix[i].getPos();
72  _I(i) = pix[i].getFlux();
73  double mask = exp(-std::norm(pix[i].getPos()/_sigma)/2.);
74  _W(i) = pix[i].getInverseSigma()*mask;
75  }
76  }
77 
78  void CrudeSolver::calculateF(const DVector& x, DVector& f) const
79  {
80  Assert(x.size() == 4);
81  Assert(f.size() == _Z.size());
82 
83  if ((x-_xInit).TMV_subVector(0,3).TMV_normInf() > 2.) {
84  f = 2.e10*_f1;
85  f.TMV_addToAll(1.);
86  return;
87  }
88  if (x(3) < 0.) {
89  f = 2.e10*_f1;
90  f.TMV_addToAll(1.);
91  return;
92  }
93 
94  std::complex<double> zc(x[0],x[1]);
95  double mu = x[2];
96  double I0 = _I1 * x[3];
97 
98  double m0 = exp(-mu)/_sigma;
99  // z' = m*(z-zc)
100  // = m*z - m*zc
101  // x' = m x - m xc
102  // y' = m y - m yc
103  _Z1 = m0*_Z;
104  _Z1.TMV_addToAll(-m0*zc);
105 
106  const int nPix = _Z.size();
107  for(int i=0;i<nPix;++i) {
108  double rsq = std::norm(_Z1[i]);
109  _Rsq[i] = rsq;
110  _E[i] = exp(-rsq/2.);
111  }
112  _f1 = _I - I0*_E;
113 #ifdef USE_TMV
114  _f1 = _W * _f1;
115 #else
116  _f1.array() *= _W.array();
117 #endif
118  f = _f1;
119  }
120 
122  const DVector& x, const DVector& f, DMatrix& df) const
123  {
124  //xdbg<<"Start J\n";
125  Assert(x.size() == 4);
126  Assert(f.size() == _Z.size());
127  Assert(df.TMV_colsize() == _Z.size());
128  Assert(df.TMV_rowsize() == 4);
129 
130  double mu = x[2];
131  double I0 = _I1 * x[3];
132  double m0 = exp(-mu)/_sigma;
133 
134  // fi = Wi * (Ii - I0 Ei)
135  // dfi/dI0 = -Wi Ei
136  // dfi/dmu = Wi I0 Ei (x1 dx1/dxc + y1 dy1/dxc)
137  // likewise for the other 4
138  //
139  // dx1/dxc = -m dy1/dxc = 0
140  // dx1/dyc = 0 dy1/dyc = -m
141  // dx1/dmu = -x1 dy1/dmu = -y1
142  //
143 
144 #ifdef USE_TMV
145  df.col(0) = -m0 * _Z1.realPart();
146  df.col(1) = -m0 * _Z1.imagPart();
147  df.col(2) = -_Rsq;
148  df.colRange(0,3) = I0 * DiagMatrixViewOf(_E) * df.colRange(0,3);
149  df.col(3) = -_I1 * _E;
150  df = _W * df;
151 #else
152  df.col(0).array() = _W.array() * _Z1.real().array();
153  df.col(0).array() = _E.array() * df.col(0).array();
154  df.col(0) *= -m0 * I0;
155  df.col(1).array() = _W.array() * _Z1.imag().array();
156  df.col(1).array() = _E.array() * df.col(1).array();
157  df.col(1) *= -m0 * I0;
158  df.col(2).array() = _W.array() * _Rsq.array();
159  df.col(2).array() = _E.array() * df.col(2).array();
160  df.col(2) *= -I0;
161  df.col(3).array() = _W.array() * _E.array();
162  df.col(3) *= -_I1;
163 #endif
164  }
165 
166  void Ellipse::crudeMeasure(const PixelList& pix, double sigma)
167  {
168  // We use as our initial estimate the exact value for a
169  // well-sampled, uniform-variance, undistorted Gaussian intensity pattern
170  // with no PSF convolution.
171  //
172  // That is, we assume the model:
173  //
174  // I(x,y) = I0 exp( -(|z-zc|^2/ (2 sigma'^2) )
175  // where zz = (z-zc)
176  // and sigma' = exp(mu) sigma
177  //
178  xdbg<<"Current centroid = "<<_cen<<std::endl;
179  xdbg<<"Current mu = "<<_mu<<std::endl;
180  xdbg<<"sigma = "<<sigma<<std::endl;
181  const int nPix = pix.size();
182 
183 #if 1
184  // With a weight of exp(-|z|^2/(2 sigma^2)),
185  // the weighted moments of this function are:
186  // Iz/I = zc / (1+exp(2mu))
187  // Irr/I = ( |zc|^2 + 2 exp(2mu) sigma^2 ) / (1+exp(2mu))
188 
189  std::complex<double> Iz = 0.;
190  double I = 0.;
191  double sig2 = sigma * exp(real(_mu));
192  for(int i=0;i<nPix;++i) {
193  double wt = exp(-std::norm((pix[i].getPos()-_cen)/sig2)/2.);
194  Iz += wt * pix[i].getFlux() * (pix[i].getPos()-_cen);
195  I += wt * pix[i].getFlux();
196  if (std::abs(pix[i].getPos()-_cen) < 2.)
197  xdbg<<pix[i].getPos()<<" "<<pix[i].getFlux()<<std::endl;
198  }
199  xdbg<<"Iz = "<<Iz<<", I = "<<I<<std::endl;
200 
201  // If I <= 0 then this isn't going to work. Just return and hope
202  // the regular measure method might do better.
203  if (!(I > 0.)) return;
204 
205  std::complex<double> zc = Iz / I;
206  // If zc is more than 2, something is probably wrong, so abort now.
207  if (!(std::abs(zc) < 2.)) return;
208 
209  xdbg<<"Initial offset to centroid = "<<zc<<std::endl;
210  if (isFixedCen()) {
211  xdbg<<"But centroid is fixed, so don't apply.\n";
212  zc = _cen;
213  } else {
214  zc += _cen;
215  }
216  xdbg<<"zc = "<<zc<<std::endl;
217 
218  double Irr = 0.;
219  double W = 0.;
220  I = 0.;
221  Iz = 0.;
222  for(int i=0;i<nPix;++i) {
223  double wt = exp(-std::norm((pix[i].getPos()-zc)/sig2)/2.);
224  Iz += wt * pix[i].getFlux() * (pix[i].getPos()-zc);
225  Irr += wt * pix[i].getFlux() * std::norm(pix[i].getPos()-zc);
226  I += wt * pix[i].getFlux();
227  W += wt;
228  }
229  xdbg<<"Iz = "<<Iz<<", Irr = "<<Irr<<", I = "<<I<<", W = "<<W<<std::endl;
230 
231  std::complex<double> zc1 = Iz/I;
232  double S = Irr/I - norm(zc1);
233  // S is now 2 exp(2mu) sigma^2 / (1 + exp(2mu))
234  xdbg<<"S = "<<S<<std::endl;
235  double exp2mu = S / 2. / (sig2*sig2);
236  xdbg<<"exp2mu/(1+exp2mu) = "<<exp2mu<<std::endl;
237  // It's actually exp(2mu) / (1+exp(2mu)) at this point
238  if (exp2mu < 0.2)
239  exp2mu = 0.25;
240  else if (exp2mu < 0.8)
241  exp2mu = 1./(1./exp2mu-1.); // Now it is really exp(2mu)
242  else
243  // The above formula is unstable, and probably inappropriate, since
244  // we probably have a failure of our model approximation.
245  // So just multiply it by 5 -- the correct factor for exp2mu = 0.8
246  exp2mu *= 5.;
247  xdbg<<"exp2mu = "<<exp2mu<<std::endl;
248  if (exp2mu <= 0.) exp2mu = 1.;
249 
250  double m = log(exp2mu)/2.;
251  xdbg<<"mu = "<<m<<std::endl;
252 
253  if (!isFixedCen()) zc += zc1 * (1.+exp2mu);
254  if (isFixedMu()) m = real(_mu);
255  else m += real(_mu);
256 
257  xdbg<<"Approx cen = "<<zc<<std::endl;
258  xdbg<<"Approx mu = "<<m<<std::endl;
259 
260  // I/W = I0 exp(2mu) / (1+exp(2mu)) * exp(-|zc1|^2/2sigma^2*(1+exp(2mu)))
261  double I0 = (I/W)*(1.+exp2mu)/exp2mu /
262  exp(-norm(zc1)*(1.+exp2mu)/(2.*sig2*sig2));
263 #else
264  std::complex<double> zc = _cen;
265  double m = real(_mu);
266 
267  double model = 0.;
268  double obs = 0.;
269  double minx = 1.e100, maxx = -1.e100, miny = 1.e100, maxy = -1.e100;
270  for(int i=0;i<nPix;++i) {
271  if (real(pix[i].getPos()) < minx) minx = real(pix[i].getPos());
272  if (real(pix[i].getPos()) > maxx) maxx = real(pix[i].getPos());
273  if (imag(pix[i].getPos()) < miny) miny = imag(pix[i].getPos());
274  if (imag(pix[i].getPos()) > maxy) maxy = imag(pix[i].getPos());
275  double wt = exp(-std::norm(pix[i].getPos()/sigma)/2.);
276  model += mask;
277  obs += pix[i].getFlux() * mask;
278  }
279  double pixarea = (maxx-minx)*(maxy-miny)/nPix;
280  xdbg<<"pixarea = "<<pixarea<<std::endl;
281  xdbg<<"obs = "<<obs<<std::endl;
282  xdbg<<"model = "<<model<<std::endl;
283  double I0 = 2.*obs / model;
284 #endif
285 
286  xdbg<<"Initial I0 estimate = "<<I0<<std::endl;
287  if (I0 < 1.e-6) {
288  xdbg<<"Warning: small or negative I0: "<<I0<<" -- Use 1.0\n";
289  I0 = 1.;
290  }
291 
292  //if (std::abs(zc) > 1.0) zc /= std::abs(zc);
293  //if (std::abs(m) > 1.0) m /= std::abs(m);
294  DVector x(4);
295  x[0] = std::real(zc); x[1] = std::imag(zc);
296  x[2] = m;
297  x[3] = 1.;
298  DVector f(nPix);
299  CrudeSolver s(pix,sigma,I0,x);
300 
301  s.useHybrid();
302  s.setTol(1.e-4,1.e-8);
303  s.setMinStep(1.e-15);
304  s.setTau(1.0);
305  s.setDelta0(0.05);
306 #ifdef __PGI
307  s.noUseCholesky();
308 #endif
309  if (XDEBUG) s.setOutput(*dbgout);
310  xdbg<<"Before CrudeSolver: x = "<<EIGEN_Transpose(x)<<std::endl;
311  s.solve(x,f);
312  xdbg<<"After CrudeSolver: x = "<<EIGEN_Transpose(x)<<std::endl;
313  s.setFTol(1.e-4 * (1.e-4 + std::abs(x[3])));
314  s.solve(x,f);
315  xdbg<<"After 2nd CrudeSolver: x = "<<EIGEN_Transpose(x)<<std::endl;
316 
317  std::complex<double> cenNew(x[0],x[1]);
318  double muNew = x[2];
319 
320  if (std::abs(cenNew-_cen) > 2.) {
321  dbg<<"Warning: large centroid shift in CrudeMeasure\n";
322  dbg<<"Old centroid = "<<_cen<<", new centroid = "<<cenNew<<std::endl;
323  cenNew = _cen + 2.*(cenNew - _cen)/std::abs(cenNew-_cen);
324  dbg<<"Scaling back to "<<cenNew<<std::endl;
325  }
326 
327  if (std::abs(muNew-real(_mu)) > 2.) {
328  dbg<<"Warning: large scale change in CrudeMeasure\n";
329  dbg<<"Old mu = "<<_mu<<", new mu = "<<muNew<<std::endl;
330  muNew = real(_mu) + 2.*(muNew - real(_mu))/std::abs(muNew-real(_mu));
331  dbg<<"Scaling back to "<<muNew<<std::endl;
332  }
333 
334  xdbg<<"Crude cen = "<<cenNew<<std::endl;
335  xdbg<<"Crude mu = "<<muNew<<std::endl;
336 
337  if (!_isFixedCen) _cen = cenNew;
338  if (!_isFixedMu) _mu = muNew;
339  }
340 
342  const std::vector<PixelList>& pix, double sigma)
343  {
344  int nPix = 0;
345  const int nPixList = pix.size();
346  for(int i=0;i<nPixList;++i) nPix += pix[i].size();
347  PixelList allPix(nPix);
348  for(int i=0,k=0;i<nPixList;++i)
349  for(size_t j=0;j<pix[i].size();++j,++k)
350  allPix[k] = pix[i][j];
351  crudeMeasure(allPix,sigma);
352  }
353 
354  void Ellipse::peakCentroid(const PixelList& pix, double maxR)
355  {
356  double IPeak = 0.;
357  std::complex<double> zPeak = 0.;
358  const int nPix = pix.size();
359  for(int i=0;i<nPix;++i) if (std::abs(pix[i].getPos()) < maxR) {
360  if (pix[i].getFlux() > IPeak) {
361  zPeak = pix[i].getPos();
362  IPeak = pix[i].getFlux();
363  }
364  }
365  _cen = zPeak;
366  }
367 
368 }}}}
const bool XDEBUG
Definition: dbg.h:63
double _sigma
Definition: SincCoeffs.cc:196
virtual void setTau(double tau)
Definition: NLSolver.h:400
virtual bool solve(DVector &x, DVector &f) const
Definition: NLSolver.cc:1620
T norm(const T &x)
Definition: Integrate.h:191
virtual void setOutput(std::ostream &os)
Definition: NLSolver.h:410
void crudeMeasure(const PixelList &pix, double sigma)
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
def log
Definition: log.py:85
void peakCentroid(const PixelList &pix, double maxr)
void calculateJ(const DVector &x, const DVector &f, DMatrix &df) const
virtual void setTol(double fTol, double gTol)
Definition: NLSolver.h:395
virtual void setDelta0(double delta0)
Definition: NLSolver.h:401
virtual void setFTol(double fTol)
Definition: NLSolver.h:393
int x
virtual void setMinStep(double minStep)
Definition: NLSolver.h:398
CrudeSolver(const PixelList &pix, double sigma, double I1, DVector &xinit)
Definition: CrudeMeasure.cc:62
tuple m
Definition: lsstimport.py:48
#define xdbg
Definition: dbg.h:70
void calculateF(const DVector &x, DVector &f) const
Definition: CrudeMeasure.cc:78
#define dbg
Definition: dbg.h:69
T real(const T &x)
Definition: Integrate.h:193
std::ostream *const dbgout
Definition: dbg.h:64
#define EIGEN_Transpose(m)
Definition: MyMatrix.h:211
#define Assert(x)
Definition: dbg.h:73