LSSTApplications  10.0+286,10.0+36,10.0+46,10.0-2-g4f67435,10.1+152,10.1+37,11.0,11.0+1,11.0-1-g47edd16,11.0-1-g60db491,11.0-1-g7418c06,11.0-2-g04d2804,11.0-2-g68503cd,11.0-2-g818369d,11.0-2-gb8b8ce7
LSSTDataManagementBasePackage
Public Member Functions | Private Member Functions | Private Attributes | List of all members
lsst::meas::algorithms::shapelet::Ellipse Class Reference

#include <Ellipse.h>

Public Member Functions

 Ellipse ()
 
 Ellipse (std::complex< double > cen, std::complex< double > gamma, std::complex< double > mu)
 
 Ellipse (double vals[])
 
 Ellipse (const Ellipse &e2)
 
Ellipse & operator= (const Ellipse &e2)
 
bool measure (const std::vector< PixelList > &pix, int order, int order2, int maxm, double sigma, long &flag, double thresh, DMatrix *cov=0)
 
bool measure (const std::vector< PixelList > &pix, const std::vector< BVec > &psf, int order, int order2, int maxm, double sigma, long &flag, double thresh, DMatrix *cov=0)
 
bool findRoundFrame (const BVec &b, bool psf, int galOrder2, double thresh, long &flag, DMatrix *cov=0)
 
void crudeMeasure (const PixelList &pix, double sigma)
 
void crudeMeasure (const std::vector< PixelList > &pix, double sigma)
 
void peakCentroid (const PixelList &pix, double maxr)
 
bool measureShapelet (const std::vector< PixelList > &pix, BVec &bret, int order, int order2, int maxm, DMatrix *bcov=0) const
 
bool measureShapelet (const std::vector< PixelList > &pix, const std::vector< BVec > &psf, BVec &bret, int order, int order2, int maxm, DMatrix *bcov=0) const
 
bool altMeasureShapelet (const std::vector< PixelList > &pix, const std::vector< BVec > &psf, BVec &bret, int order, int order2, double pixScale, DMatrix *bcov=0) const
 
bool altMeasureShapelet (const std::vector< PixelList > &pix, BVec &bret, int order, int order2, double pixScale, DMatrix *bcov=0) const
 
void preShiftBy (const std::complex< double > &cen, const std::complex< double > &gamma, const std::complex< double > &mu)
 
void postShiftBy (const std::complex< double > &cen, const std::complex< double > &gamma, const std::complex< double > &mu)
 
void removeRotation ()
 
std::complex< double > getCen () const
 
std::complex< double > getGamma () const
 
double getMu () const
 
double getTheta () const
 
void setCen (const std::complex< double > &cen)
 
void setGamma (const std::complex< double > &gamma)
 
void setMu (double mu)
 
void setTheta (double theta)
 
void fixCen ()
 
void fixGam ()
 
void fixMu ()
 
void unfixCen ()
 
void unfixGam ()
 
void unfixMu ()
 
bool isFixedCen () const
 
bool isFixedGamma () const
 
bool isFixedMu () const
 
void doTimings ()
 
void write (std::ostream &os) const
 

Private Member Functions

bool doMeasure (const std::vector< PixelList > &pix, const std::vector< BVec > *psf, int order, int order2, int maxm, double sigma, long &flag, double thresh, DMatrix *cov=0)
 
bool doMeasureShapelet (const std::vector< PixelList > &pix, const std::vector< BVec > *psf, BVec &bret, int order, int order2, int maxm, DMatrix *bcov=0) const
 
bool doAltMeasureShapelet (const std::vector< PixelList > &pix, const std::vector< BVec > *psf, BVec &bret, int order, int order2, double pixScale, DMatrix *bcov=0) const
 

Private Attributes

std::complex< double > _cen
 
std::complex< double > _gamma
 
std::complex< double > _mu
 
bool _isFixedCen
 
bool _isFixedGamma
 
bool _isFixedMu
 
bool _shouldDoTimings
 

Detailed Description

Definition at line 16 of file Ellipse.h.

Constructor & Destructor Documentation

lsst::meas::algorithms::shapelet::Ellipse::Ellipse ( )
inline

Definition at line 21 of file Ellipse.h.

21  :
22  _cen(0.), _gamma(0.), _mu(0.),
23  _isFixedCen(false), _isFixedGamma(false), _isFixedMu(false),
24  _shouldDoTimings(false) {}
lsst::meas::algorithms::shapelet::Ellipse::Ellipse ( std::complex< double >  cen,
std::complex< double >  gamma,
std::complex< double >  mu 
)
inline

Definition at line 26 of file Ellipse.h.

27  :
28  _cen(cen), _gamma(gamma), _mu(mu),
29  _isFixedCen(false), _isFixedGamma(false), _isFixedMu(false),
30  _shouldDoTimings(false) {}
lsst::meas::algorithms::shapelet::Ellipse::Ellipse ( double  vals[])
inline

Definition at line 32 of file Ellipse.h.

32  :
33  _cen(vals[0],vals[1]), _gamma(vals[2],vals[3]), _mu(vals[4]),
34  _isFixedCen(false), _isFixedGamma(false), _isFixedMu(false),
35  _shouldDoTimings(false) {}
lsst::meas::algorithms::shapelet::Ellipse::Ellipse ( const Ellipse &  e2)
inline

Definition at line 39 of file Ellipse.h.

39  :
40  _cen(e2.getCen()), _gamma(e2.getGamma()),
41  _mu(e2.getMu(),e2.getTheta()),
42  _isFixedCen(false), _isFixedGamma(false), _isFixedMu(false),
43  _shouldDoTimings(false) {}

Member Function Documentation

bool lsst::meas::algorithms::shapelet::Ellipse::altMeasureShapelet ( const std::vector< PixelList > &  pix,
const std::vector< BVec > &  psf,
BVec &  bret,
int  order,
int  order2,
double  pixScale,
DMatrix *  bcov = 0 
) const

Definition at line 619 of file Ellipse.cc.

623  { return doAltMeasureShapelet(pix,&psf,b,order,order2,pixScale,bCov); }
lsst::meas::algorithms::shapelet::Ellipse::doAltMeasureShapelet
bool doAltMeasureShapelet(const std::vector< PixelList > &pix, const std::vector< BVec > *psf, BVec &bret, int order, int order2, double pixScale, DMatrix *bcov=0) const
Definition: Ellipse.cc:445
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
bool lsst::meas::algorithms::shapelet::Ellipse::altMeasureShapelet ( const std::vector< PixelList > &  pix,
BVec &  bret,
int  order,
int  order2,
double  pixScale,
DMatrix *  bcov = 0 
) const

Definition at line 625 of file Ellipse.cc.

628  { return doAltMeasureShapelet(pix,0,b,order,order2,pixScale,bCov); }
lsst::meas::algorithms::shapelet::Ellipse::doAltMeasureShapelet
bool doAltMeasureShapelet(const std::vector< PixelList > &pix, const std::vector< BVec > *psf, BVec &bret, int order, int order2, double pixScale, DMatrix *bcov=0) const
Definition: Ellipse.cc:445
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
void lsst::meas::algorithms::shapelet::Ellipse::crudeMeasure ( const PixelList &  pix,
double  sigma 
)

Definition at line 166 of file CrudeMeasure.cc.

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  }
XDEBUG
const bool XDEBUG
Definition: dbg.h:63
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
lsst.log.log.log
def log
Definition: log.py:85
x
double x
Definition: ChebyshevBoundedField.cc:305
lsstimport.m
tuple m
Definition: lsstimport.py:48
xdbg
#define xdbg
Definition: dbg.h:70
dbg
#define dbg
Definition: dbg.h:69
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
dbgout
std::ostream *const dbgout
Definition: dbg.h:64
EIGEN_Transpose
#define EIGEN_Transpose(m)
Definition: MyMatrix.h:211
void lsst::meas::algorithms::shapelet::Ellipse::crudeMeasure ( const std::vector< PixelList > &  pix,
double  sigma 
)

Definition at line 341 of file CrudeMeasure.cc.

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  }
lsst::meas::algorithms::shapelet::Ellipse::crudeMeasure
void crudeMeasure(const PixelList &pix, double sigma)
Definition: CrudeMeasure.cc:166
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
bool lsst::meas::algorithms::shapelet::Ellipse::doAltMeasureShapelet ( const std::vector< PixelList > &  pix,
const std::vector< BVec > *  psf,
BVec &  bret,
int  order,
int  order2,
double  pixScale,
DMatrix *  bcov = 0 
) const
private

Definition at line 445 of file Ellipse.cc.

449  {
450  xdbg<<"Start AltMeasureShapelet: order = "<<order<<std::endl;
451  xdbg<<"b.order, sigma = "<<b.getOrder()<<", "<<b.getSigma()<<std::endl;
452  xdbg<<"el = "<<*this<<std::endl;
453  xdbg<<"pixScale = "<<pixScale<<std::endl;
454  const int nExp = pix.size();
455 
456  const int bsize = (order+1)*(order+2)/2;
457  const int bsize2 = (order2+1)*(order2+2)/2;
458  b.vec().setZero();
459  double sigma = b.getSigma();
460 
461  double gsq = std::norm(_gamma);
462  Assert(gsq < 1.);
463  std::complex<double> m = exp(-_mu)/sqrt(1.-gsq);
464 
465  std::auto_ptr<DMatrix> cov1;
466  if (bCov) {
467  bCov->setZero();
468  cov1.reset(new DMatrix(bsize2,bsize2));
469  }
470 
471  for(int k=0;k<nExp;++k) {
472  const int nPix = pix[k].size();
473  DVector I(nPix);
474  CDVector Z(nPix);
475 
476  double sigma_obs =
477  psf ?
478  sqrt(pow(sigma,2)+pow((*psf)[k].getSigma(),2)) :
479  sigma;
480  dbg<<"sigma_obs = "<<sigma_obs<<std::endl;
481 
482  DDiagMatrix normD(bsize2);
483  double val = pow(pixScale/sigma_obs,2) * exp(-2.*real(_mu)) / nExp;
484  for(int n=0,i=0;n<=order2;++n) {
485  for(int p=n,q=0;p>=q;--p,++q,++i) {
486  if (p!=q) { normD(i) = val/2.; normD(++i) = val/2.; }
487  else normD(i) = val;
488  }
489  }
490 
491  for(int i=0;i<nPix;++i) {
492  I(i) = pix[k][i].getFlux();
493  std::complex<double> z1 = pix[k][i].getPos();
494  std::complex<double> z2 = m*((z1-_cen) - _gamma*conj(z1-_cen));
495  Z(i) = z2 / sigma_obs;
496  }
497 
498  DMatrix A(nPix,bsize2);
499  makePsi(A,TMV_vview(Z),order2);
500 
501  // b = C^-1 normD At I
502  DVector b1 = A.transpose() * I;
503  dbg<<"b1 = "<<b1<<std::endl;
504  b1 = normD EIGEN_asDiag() * b1;
505  dbg<<"b1 => "<<b1<<std::endl;
506 
507  if (bCov) {
508  // cov(I) = diag(sigma_i^2)
509  // Propagate this through to cov(b1)
510  // b1 = normD At I
511  // cov1 = normD At cov(I) A normD
512  DDiagMatrix W(nPix);
513  for(int i=0;i<nPix;++i) {
514  W(i) = 1./pix[k][i].getInverseSigma();
515  }
516  DMatrix temp = normD * A.transpose() * W;
517  *cov1 = temp * temp.transpose();
518  }
519 
520  if (psf) {
521  xdbg<<"psf = "<<(*psf)[k]<<std::endl;
522  int psforder = (*psf)[k].getOrder();
523  int newpsforder = std::max(psforder,order2);
524  xdbg<<"psforder = "<<psforder<<std::endl;
525  xdbg<<"newpsforder = "<<newpsforder<<std::endl;
526  const double psfsigma = (*psf)[k].getSigma();
527  xdbg<<"sigma = "<<psfsigma<<std::endl;
528  BVec newpsf(newpsforder,psfsigma);
529  int psfsize = (*psf)[k].size();
530  int newpsfsize = newpsf.size();
531  xdbg<<"psfsize = "<<psfsize<<std::endl;
532  bool setnew = false;
533  if (_gamma != 0.) {
534  DMatrix S(newpsfsize,psfsize);
535  calculateGTransform(_gamma,newpsforder,psforder,S);
536  newpsf.vec() = S * (*psf)[k].vec();
537  setnew = true;
538  xdbg<<"newpsf = "<<newpsf<<std::endl;
539  }
540  if (real(_mu) != 0.) {
541  if (setnew) {
542  DMatrix D(newpsfsize,newpsfsize);
543  calculateMuTransform(real(_mu),newpsforder,D);
544  newpsf.vec() = D * newpsf.vec();
545  } else {
546  DMatrix D(newpsfsize,psfsize);
547  calculateMuTransform(real(_mu),newpsforder,psforder,D);
548  newpsf.vec() = D * (*psf)[k].vec();
549  setnew = true;
550  }
551  newpsf.vec() *= exp(2.*real(_mu));
552  xdbg<<"newpsf => "<<newpsf<<std::endl;
553  }
554  if (imag(_mu) != 0.) {
555  if (setnew) {
556 #ifdef USE_TMV
557  DBandMatrix R(newpsfsize,newpsfsize,1,1);
558 #else
559  DBandMatrix R(newpsfsize,newpsfsize);
560 #endif
561  calculateThetaTransform(imag(_mu),newpsforder,R);
562  newpsf.vec() = R * newpsf.vec();
563  } else {
564 #ifdef USE_TMV
565  DBandMatrix R(newpsfsize,psfsize,1,1);
566 #else
567  DBandMatrix R(newpsfsize,psfsize);
568 #endif
569  calculateThetaTransform(imag(_mu),newpsforder,psforder,R);
570  newpsf.vec() = R * (*psf)[k].vec();
571  setnew = true;
572  }
573  xdbg<<"newpsf => "<<newpsf<<std::endl;
574  }
575 
576  DMatrix C(bsize2,bsize2);
577  if (setnew) {
578  calculatePsfConvolve(newpsf,order2,b.getSigma(),C);
579  } else {
580  calculatePsfConvolve((*psf)[k],order2,b.getSigma(),C);
581  }
582 #ifdef USE_TMV
583  b1 /= C;
584 #else
585  b1 = C.lu().solve(b1);
586 #endif
587  dbg<<"b1 /= C => "<<b1<<std::endl;
588  if (bCov) *cov1 = C.inverse() * (*cov1) * C;
589  }
590  b.vec() += b1.TMV_subVector(0,bsize);
591  if (bCov) *bCov += cov1->TMV_subMatrix(0,bsize,0,bsize);
592  }
593  if (!(b(0) > 0.)) {
594  dbg<<"Calculated b vector has negative b(0):\n";
595  dbg<<"b = "<<b<<std::endl;
596  return false;
597  }
598  if (order < b.getOrder()) {
599  xdbg<<"Need to zero the rest of b\n";
600  xdbg<<"bsize = "<<bsize<<" b.size = "<<b.size()<<std::endl;
601  // Zero out the rest of the shapelet vector:
602  b.vec().TMV_subVector(bsize,b.size()).setZero();
603  }
604  xdbg<<"Done (alt) measure Shapelet\n";
605  return true;
606  }
TMV_vview
#define TMV_vview(v)
Definition: MyMatrix.h:202
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
lsst::meas::algorithms::shapelet::calculateGTransform
void calculateGTransform(std::complex< double > g, int order1, int order2, DMatrix &S)
Definition: BVec.cc:468
val
bool val
Definition: DevicePreference.cc:55
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
lsst::meas::algorithms::shapelet::makePsi
void makePsi(DMatrix &psi, CDVectorView z, int order, const DVector *coeff=0)
Definition: PsiHelper.cc:226
lsst::meas::algorithms::shapelet::calculateThetaTransform
void calculateThetaTransform(double theta, int order1, int order2, DBandMatrix &R)
Definition: BVec.cc:419
EIGEN_asDiag
#define EIGEN_asDiag()
Definition: MyMatrix.h:180
lsstimport.m
tuple m
Definition: lsstimport.py:48
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
xdbg
#define xdbg
Definition: dbg.h:70
lsst::meas::algorithms::shapelet::calculatePsfConvolve
void calculatePsfConvolve(const BVec &bpsf, int order1, int order2, double sigma, DMatrix &C)
Definition: BVec.cc:675
dbg
#define dbg
Definition: dbg.h:69
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
Assert
#define Assert(x)
Definition: dbg.h:73
lsst::meas::algorithms::shapelet::calculateMuTransform
void calculateMuTransform(double mu, int order1, int order2, DMatrix &D)
Definition: BVec.cc:227
bool lsst::meas::algorithms::shapelet::Ellipse::doMeasure ( const std::vector< PixelList > &  pix,
const std::vector< BVec > *  psf,
int  order,
int  order2,
int  maxm,
double  sigma,
long &  flag,
double  thresh,
DMatrix *  cov = 0 
)
private

Definition at line 39 of file Ellipse_meas.cc.

44  {
45  dbg<<"Start DoMeasure: galOrder = "<<galOrder<<", psf = "<<bool(psf)<<std::endl;
46  dbg<<"fix = "<<_isFixedCen<<" "<<_isFixedGamma<<" "<<_isFixedMu<<std::endl;
47  dbg<<"Thresh = "<<thresh<<std::endl;
48  int npix = 0;
49  for(size_t i=0;i<pix.size();++i) {
50  xdbg<<"npix["<<i<<"] = "<<pix[i].size()<<std::endl;
51  npix += pix[i].size();
52  }
53 
54  int galSize = (galOrder+1)*(galOrder+2)/2;
55  if (npix <= galSize) {
56  dbg<<"Too few pixels ("<<npix<<") for given gal_order. \n";
57  return false;
58  }
59 
60  BVec b(galOrder,sigma);
61 
62  if (!doMeasureShapelet(pix,psf,b,galOrder,galOrder2,maxm)) {
63  xdbg<<"Could not measure a shapelet vector.\n";
64  return false;
65  }
66  if (b(0) <= 0) {
67  xdbg<<"Bad flux in measured shapelet\n";
68  return false;
69  }
70  xdbg<<"b = "<<b<<std::endl;
71 
72  return findRoundFrame(b,psf,galOrder2,thresh,flag,cov);
73  }
lsst::meas::algorithms::shapelet::Ellipse::findRoundFrame
bool findRoundFrame(const BVec &b, bool psf, int galOrder2, double thresh, long &flag, DMatrix *cov=0)
Definition: Ellipse_meas.cc:75
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
lsst::meas::algorithms::shapelet::Ellipse::doMeasureShapelet
bool doMeasureShapelet(const std::vector< PixelList > &pix, const std::vector< BVec > *psf, BVec &bret, int order, int order2, int maxm, DMatrix *bcov=0) const
Definition: Ellipse.cc:186
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
xdbg
#define xdbg
Definition: dbg.h:70
dbg
#define dbg
Definition: dbg.h:69
bool lsst::meas::algorithms::shapelet::Ellipse::doMeasureShapelet ( const std::vector< PixelList > &  pix,
const std::vector< BVec > *  psf,
BVec &  bret,
int  order,
int  order2,
int  maxm,
DMatrix *  bcov = 0 
) const
private

Definition at line 186 of file Ellipse.cc.

190  {
191  xdbg<<"Start MeasureShapelet: order = "<<order<<std::endl;
192  xdbg<<"b.order, sigma = "<<b.getOrder()<<", "<<b.getSigma()<<std::endl;
193  xdbg<<"el = "<<*this<<std::endl;
194  if (maxm < 0 || maxm > order) maxm = order;
195  xdbg<<"order = "<<order<<','<<order2<<','<<maxm<<std::endl;
196 
197  // ( u )' = exp(-mu)/sqrt(1-gsq) ( 1-g1 -g2 ) ( u-uc )
198  // ( v ) ( -g2 1+g1 ) ( v-vc )
199  //
200  // z' = u' + I v'
201  // = exp(-mu)/sqrt(1-gsq)
202  // [ (1-g1)(u-uc)-g2(v-vc)-Ig2(u-uc)+I(1+g1)(v-vc) ]
203  // = exp(-mu)/sqrt(1-gsq) [ z-zc - g1(z-zc)* -Ig2(z-zc)* ]
204  // = exp(-mu)/sqrt(1-gsq) ( z-zc - g (z-zc)* )
205 
206  double sigma = b.getSigma();
207  double gsq = std::norm(_gamma);
208  Assert(gsq < 1.);
209 
210  std::complex<double> m = exp(-_mu)/sqrt(1.-gsq);
211 
212  int nTot = 0;
213  const int nExp = pix.size();
214  for(int i=0;i<nExp;++i) nTot += pix[i].size();
215  dbg<<"ntot = "<<nTot<<" in "<<nExp<<" images\n";
216 
217  DVector I(nTot);
218  DVector W(nTot);
219  CDVector Z(nTot);
220 
221  for(int k=0,n=0;k<nExp;++k) {
222  double sigma_obs =
223  psf ?
224  sqrt(pow(sigma,2)+pow((*psf)[k].getSigma(),2)) :
225  sigma;
226  dbg<<"sigma_obs["<<k<<"] = "<<sigma_obs<<std::endl;
227 
228  const int nPix = pix[k].size();
229  for(int i=0;i<nPix;++i,++n) {
230  I(n) = pix[k][i].getFlux()*pix[k][i].getInverseSigma();
231  W(n) = pix[k][i].getInverseSigma();
232  std::complex<double> z1 = pix[k][i].getPos();
233  std::complex<double> z2 = m*((z1-_cen) - _gamma*conj(z1-_cen));
234  Z(n) = z2 / sigma_obs;
235  }
236  }
237  //xdbg<<"Z = "<<Z<<std::endl;
238  //xdbg<<"I = "<<I<<std::endl;
239  //xdbg<<"W = "<<W<<std::endl;
240 
241  int bsize = (order+1)*(order+2)/2;
242  xdbg<<"bsize = "<<bsize<<std::endl;
243  int bsize2 = (order2+1)*(order2+2)/2;
244  xdbg<<"bsize2 = "<<bsize2<<std::endl;
245 
246 #ifdef USE_TMV
247  // Figure out a permutation that puts all the m <= maxm first.
248  // I don't know how to do this with Eigen, but I don't really
249  // use maxm < order on a regular basis, so only implement this for TMV.
250  tmv::Permutation P(bsize);
251  int msize = bsize;
252  if (maxm < order) {
253  tmv::Vector<double> mvals(bsize);
254  for(int n=0,k=0;n<=order;++n) {
255  for(int m=n;m>=0;m-=2) {
256  mvals[k++] = m;
257  if (m > 0) mvals[k++] = m;
258  }
259  }
260  dbg<<"mvals = "<<mvals<<std::endl;
261  mvals.sort(P);
262  dbg<<"mvals => "<<mvals<<std::endl;
263  // 00 10 10 20 20 11 30 30 21 21 40 40 31 31 22 50 50 41 41 32 32
264  // n = 0 1 2 3 4 5 6
265  // 0size = 1 1 2 2 3 3 4 = (n)/2 + 1
266  // 1size = 1 3 4 6 7 9 10 = (3*(n-1))/2 + 3
267  // 2size = 1 3 6 8 11 13 16 = (5*(n-2))/2 + 6
268  // 3size = 1 3 6 10 13 17 20 = (7*(n-3))/2 + 10
269  // nsize = (n+1)*(n+2)/2
270  // msize = (m+1)*(m+2)/2 + (2*m+1)*(n-m)/2
271  msize = (maxm+1)*(maxm+2)/2 + (2*maxm+1)*(order-maxm)/2;
272  dbg<<"msize = "<<msize<<std::endl;
273  dbg<<"mvals["<<msize-1<<"] = "<<mvals[msize-1]<<std::endl;
274  dbg<<"mvals["<<msize<<"] = "<<mvals[msize]<<std::endl;
275  Assert(mvals[msize-1] == maxm);
276  Assert(mvals[msize] == maxm+1);
277  }
278 #endif
279 
280  DMatrix A(nTot,bsize);
281  Assert(nTot >= bsize); // Should have been addressed by calling routine.
282  //xdbg<<"A = "<<A<<std::endl;
283 
284  if (psf) {
285  for(int k=0,n=0,nx;k<nExp;++k,n=nx) {
286  xdbg<<"psf = "<<(*psf)[k]<<std::endl;
287  int psforder = (*psf)[k].getOrder();
288  int newpsforder = std::max(psforder,order2);
289  xdbg<<"psforder = "<<psforder<<std::endl;
290  xdbg<<"newpsforder = "<<newpsforder<<std::endl;
291  const double psfsigma = (*psf)[k].getSigma();
292  xdbg<<"sigma = "<<psfsigma<<std::endl;
293  BVec newpsf(newpsforder,psfsigma);
294  int psfsize = (*psf)[k].size();
295  int newpsfsize = newpsf.size();
296  xdbg<<"psfsize = "<<psfsize<<std::endl;
297  bool setnew = false;
298  if (_gamma != 0.) {
299  DMatrix S(newpsfsize,psfsize);
300  calculateGTransform(_gamma,newpsforder,psforder,S);
301  newpsf.vec() = S * (*psf)[k].vec();
302  setnew = true;
303  xdbg<<"newpsf = "<<newpsf<<std::endl;
304  }
305  if (real(_mu) != 0.) {
306  if (setnew) {
307  DMatrix D(newpsfsize,newpsfsize);
308  calculateMuTransform(real(_mu),newpsforder,D);
309  newpsf.vec() = D * newpsf.vec();
310  } else {
311  DMatrix D(newpsfsize,psfsize);
312  calculateMuTransform(real(_mu),newpsforder,psforder,D);
313  newpsf.vec() = D * (*psf)[k].vec();
314  setnew = true;
315  }
316  newpsf.vec() *= exp(2.*real(_mu));
317  xdbg<<"newpsf => "<<newpsf<<std::endl;
318  }
319  if (imag(_mu) != 0.) {
320  if (setnew) {
321 #ifdef USE_TMV
322  DBandMatrix R(newpsfsize,newpsfsize,1,1);
323 #else
324  DBandMatrix R(newpsfsize,newpsfsize);
325 #endif
326  calculateThetaTransform(imag(_mu),newpsforder,R);
327  newpsf.vec() = R * newpsf.vec();
328  } else {
329 #ifdef USE_TMV
330  DBandMatrix R(newpsfsize,psfsize,1,1);
331 #else
332  DBandMatrix R(newpsfsize,psfsize);
333 #endif
334  calculateThetaTransform(imag(_mu),newpsforder,psforder,R);
335  newpsf.vec() = R * (*psf)[k].vec();
336  setnew = true;
337  }
338  xdbg<<"newpsf => "<<newpsf<<std::endl;
339  }
340  DMatrix C(bsize2,bsize);
341  if (setnew) {
342  calculatePsfConvolve(newpsf,order2,order,b.getSigma(),C);
343  } else {
344  calculatePsfConvolve((*psf)[k],order2,order,b.getSigma(),C);
345  }
346 
347  const int nPix = pix[k].size();
348  nx = n+nPix;
349  DMatrix A1(nPix,bsize2);
350  makePsi(A1,Z.TMV_subVector(n,nx),order2,&W);
351  TMV_rowRange(A,n,nx) = A1 * C;
352  }
353  } else {
354  makePsi(A,TMV_vview(Z),order,&W);
355  }
356  const double MAX_CONDITION = 1.e8;
357 
358 #ifdef USE_TMV
359  A *= P.transpose();
360  tmv::MatrixView<double> Am = A.colRange(0,msize);
361  // I used to use SVD for this, rather than the QRP decomposition.
362  // But QRP is significantly faster, and the condition estimate
363  // produced is good enough for what we need here.
364  // TODO: I need to add a real condition estimator to division methods
365  // other than SVD in TMV. LAPACK has this functionality...
366  Am.saveDiv();
367  Am.divideUsing(tmv::QRP);
368  Am.qrpd();
369  xdbg<<"R diag = "<<Am.qrpd().getR().diag()<<std::endl;
370  if (
371 #if TMV_VERSION_AT_LEAST(0,65)
372  Am.qrpd().isSingular() ||
373 #else
374  // Fixed a bug in the isSingular function in v0.65.
375  // Until that is the standard TMV version for DES, use this instead:
376  (Am.qrpd().getR().minAbs2Element() <
377  1.e-16 * Am.qrpd().getR().maxAbs2Element()) ||
378 #endif
379  DiagMatrixViewOf(Am.qrpd().getR().diag()).condition() >
380  MAX_CONDITION) {
381  dbg<<"Poor condition in MeasureShapelet: \n";
382  dbg<<"R diag = "<<Am.qrpd().getR().diag()<<std::endl;
383  dbg<<"condition = "<<
384  DiagMatrixViewOf(Am.qrpd().getR().diag()).condition()<<
385  std::endl;
386  dbg<<"det = "<<Am.qrpd().det()<<std::endl;
387  return false;
388  }
389  b.vec().subVector(0,msize) = I/Am;
390  xdbg<<"b = "<<b<<std::endl;
391  xdbg<<"Norm(I) = "<<Norm(I)<<std::endl;
392  xdbg<<"Norm(Am*b) = "<<Norm(Am*b.vec().subVector(0,msize))<<std::endl;
393  xdbg<<"Norm(I-Am*b) = "<<Norm(I-Am*b.vec().subVector(0,msize))<<std::endl;
394  b.vec().subVector(msize,bsize).setZero();
395  b.vec() = P.transpose() * b.vec();
396  xdbg<<"b => "<<b<<std::endl;
397 
398  if (bCov) {
399  bCov->setZero();
400  Am.makeInverseATA(bCov->subMatrix(0,msize,0,msize));
401  *bCov = P.transpose() * (*bCov) * P;
402  }
403 #else
404  //const double sqrtEps = sqrt(std::numeric_limits<double>::epsilon());
405  Eigen::JacobiSVD<DMatrix> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
406  const DVector& svd_s = svd.singularValues();
407  double max = svd_s(0);
408  double min = svd_s(svd_s.size()-1);
409  if (max > MAX_CONDITION * min) {
410  xdbg<<"Poor condition in MeasureShapelet: \n";
411  xdbg<<"svd = "<<svd_s.transpose()<<std::endl;
412  xdbg<<"condition = "<<max/min<<std::endl;
413  return false;
414  }
415  // USVtb = I
416  // b = VS^-1UtI
417  DVector temp = svd.matrixU().transpose() * I;
418  temp = svd_s.array().inverse().matrix().asDiagonal() * temp;
419  b.vec().TMV_subVector(0,bsize) = svd.matrixV() * temp;
420 
421  if (bCov) {
422  bCov->setZero();
423  // (AtA)^-1 = (VSUt USVt)^-1 = (V S^2 Vt)^-1 = V S^-2 Vt
424  DMatrix temp2 =
425  svd_s.array().square().inverse().matrix().asDiagonal() * svd.matrixV().transpose();
426  bCov->TMV_subMatrix(0,bsize,0,bsize) = svd.matrixV() * temp2;
427  }
428 #endif
429  if (!(b(0) > 0.)) {
430  dbg<<"Calculated b vector has negative b(0):\n";
431  dbg<<"b = "<<b<<std::endl;
432  return false;
433  }
434  if (order < b.getOrder()) {
435  xdbg<<"Need to zero the rest of b\n";
436  xdbg<<"bsize = "<<bsize<<" b.size = "<<b.size()<<std::endl;
437  // Zero out the rest of the shapelet vector:
438  b.vec().TMV_subVector(bsize,b.size()).setZero();
439  }
440  xdbg<<"Done measure Shapelet\n";
441  xdbg<<"b = "<<b.vec()<<std::endl;
442  return true;
443  }
TMV_vview
#define TMV_vview(v)
Definition: MyMatrix.h:202
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
lsst::meas::algorithms::shapelet::calculateGTransform
void calculateGTransform(std::complex< double > g, int order1, int order2, DMatrix &S)
Definition: BVec.cc:468
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
lsst::meas::algorithms::shapelet::makePsi
void makePsi(DMatrix &psi, CDVectorView z, int order, const DVector *coeff=0)
Definition: PsiHelper.cc:226
lsst::meas::algorithms::shapelet::calculateThetaTransform
void calculateThetaTransform(double theta, int order1, int order2, DBandMatrix &R)
Definition: BVec.cc:419
lsstimport.m
tuple m
Definition: lsstimport.py:48
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
TMV_rowRange
#define TMV_rowRange(m, i1, i2)
Definition: MyMatrix.h:204
xdbg
#define xdbg
Definition: dbg.h:70
lsst::meas::algorithms::shapelet::calculatePsfConvolve
void calculatePsfConvolve(const BVec &bpsf, int order1, int order2, double sigma, DMatrix &C)
Definition: BVec.cc:675
dbg
#define dbg
Definition: dbg.h:69
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
Assert
#define Assert(x)
Definition: dbg.h:73
lsst::meas::algorithms::shapelet::calculateMuTransform
void calculateMuTransform(double mu, int order1, int order2, DMatrix &D)
Definition: BVec.cc:227
void lsst::meas::algorithms::shapelet::Ellipse::doTimings ( )
inline

Definition at line 154 of file Ellipse.h.

154 { _shouldDoTimings = true; }
bool lsst::meas::algorithms::shapelet::Ellipse::findRoundFrame ( const BVec &  b,
bool  psf,
int  galOrder2,
double  thresh,
long &  flag,
DMatrix *  cov = 0 
)

Definition at line 75 of file Ellipse_meas.cc.

78  {
79  DVector x(5);
80  DVector f(5);
81 
82  x.setZero();
83 
84  EllipseSolver3 solver(b,galOrder2,_isFixedCen,_isFixedGamma,_isFixedMu);
85 
86 #ifdef NOTHROW
87  solver.noUseCholesky();
88 #endif
89  double ftol = thresh*thresh;
90  double gtol = thresh*ftol;
91  solver.setTol(ftol,gtol);
92  solver.setMinStep(gtol*thresh);
93  solver.setOutput(*dbgout);
94  if (XDEBUG) solver.useVerboseOutput();
95  solver.setMinStep(1.e-6*gtol);
96  solver.setDelta0(0.01);
97  solver.setMaxIter(200);
98  if (psf && !_isFixedMu) {
99  solver.setDelta0(0.01);
100  solver.useDogleg();
101  solver.dontZeroB11();
102  solver.useSVD();
103  } else {
104  solver.useHybrid();
105  }
106  if (solver.solve(x,f)) {
107  dbg<<"Found good round frame:\n";
108  dbg<<"x = "<<EIGEN_Transpose(x)<<std::endl;
109  dbg<<"f = "<<EIGEN_Transpose(f)<<std::endl;
110  double f_normInf = f.TMV_normInf();
111  if (psf && !_isFixedMu && !(f_normInf < solver.getFTol())) {
112  xdbg<<"Oops, Local minimum, not real solution.\n";
113  xdbg<<"f.norm = "<<f.norm()<<std::endl;
114  xdbg<<"f.normInf = "<<f_normInf<<std::endl;
115  xdbg<<"ftol = "<<solver.getFTol()<<std::endl;
116  dbg<<"FLAG SHEAR_LOCAL_MIN\n";
117  flag |= SHEAR_LOCAL_MIN;
118  return false;
119  }
120  } else {
121  dbg<<"findRoundFrame solver failed:\n";
122  dbg<<"x = "<<EIGEN_Transpose(x)<<std::endl;
123  dbg<<"f = "<<EIGEN_Transpose(f)<<std::endl;
124  //if (XDEBUG) if (!solver.testJ(x,f,dbgout,1.e-5)) exit(1);
125  dbg<<"FLAG SHEAR_DIDNT_CONVERGE\n";
126  flag |= SHEAR_DIDNT_CONVERGE;
127  return false;
128  }
129 
130  double sigma = b.getSigma();
131  preShiftBy(std::complex<double>(x(0),x(1))*sigma,
132  std::complex<double>(x(2),x(3)),
133  x(4));
134 
135  dbg<<"ell => "<<*this<<std::endl;
136 
137  if (cov) {
138  solver.useSVD();
139  solver.getCovariance(*cov);
140  xdbg<<"cov = "<<*cov<<std::endl;
141  }
142 
143  return true;
144  }
XDEBUG
const bool XDEBUG
Definition: dbg.h:63
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
lsst::meas::algorithms::shapelet::Ellipse::preShiftBy
void preShiftBy(const std::complex< double > &cen, const std::complex< double > &gamma, const std::complex< double > &mu)
Definition: Ellipse.cc:124
x
double x
Definition: ChebyshevBoundedField.cc:305
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
xdbg
#define xdbg
Definition: dbg.h:70
dbg
#define dbg
Definition: dbg.h:69
dbgout
std::ostream *const dbgout
Definition: dbg.h:64
EIGEN_Transpose
#define EIGEN_Transpose(m)
Definition: MyMatrix.h:211
void lsst::meas::algorithms::shapelet::Ellipse::fixCen ( )
inline

Definition at line 142 of file Ellipse.h.

142 { _isFixedCen = true; }
void lsst::meas::algorithms::shapelet::Ellipse::fixGam ( )
inline

Definition at line 143 of file Ellipse.h.

143 { _isFixedGamma = true; }
void lsst::meas::algorithms::shapelet::Ellipse::fixMu ( )
inline

Definition at line 144 of file Ellipse.h.

144 { _isFixedMu = true; }
std::complex<double> lsst::meas::algorithms::shapelet::Ellipse::getCen ( ) const
inline

Definition at line 132 of file Ellipse.h.

132 { return _cen; }
std::complex<double> lsst::meas::algorithms::shapelet::Ellipse::getGamma ( ) const
inline

Definition at line 133 of file Ellipse.h.

133 { return _gamma; }
double lsst::meas::algorithms::shapelet::Ellipse::getMu ( ) const
inline

Definition at line 134 of file Ellipse.h.

134 { return real(_mu); }
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
double lsst::meas::algorithms::shapelet::Ellipse::getTheta ( ) const
inline

Definition at line 135 of file Ellipse.h.

135 { return imag(_mu); }
bool lsst::meas::algorithms::shapelet::Ellipse::isFixedCen ( ) const
inline

Definition at line 150 of file Ellipse.h.

150 { return _isFixedCen; }
bool lsst::meas::algorithms::shapelet::Ellipse::isFixedGamma ( ) const
inline

Definition at line 151 of file Ellipse.h.

151 { return _isFixedGamma; }
bool lsst::meas::algorithms::shapelet::Ellipse::isFixedMu ( ) const
inline

Definition at line 152 of file Ellipse.h.

152 { return _isFixedMu; }
bool lsst::meas::algorithms::shapelet::Ellipse::measure ( const std::vector< PixelList > &  pix,
int  order,
int  order2,
int  maxm,
double  sigma,
long &  flag,
double  thresh,
DMatrix *  cov = 0 
)

Definition at line 180 of file Ellipse.cc.

184  { return doMeasure(pix,0,order,order2,maxm,sigma,flag,thresh,cov); }
lsst::meas::algorithms::shapelet::Ellipse::doMeasure
bool doMeasure(const std::vector< PixelList > &pix, const std::vector< BVec > *psf, int order, int order2, int maxm, double sigma, long &flag, double thresh, DMatrix *cov=0)
Definition: Ellipse_meas.cc:39
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
bool lsst::meas::algorithms::shapelet::Ellipse::measure ( const std::vector< PixelList > &  pix,
const std::vector< BVec > &  psf,
int  order,
int  order2,
int  maxm,
double  sigma,
long &  flag,
double  thresh,
DMatrix *  cov = 0 
)

Definition at line 173 of file Ellipse.cc.

178  { return doMeasure(pix,&psf,order,order2,maxm,sigma,flag,thresh,cov); }
lsst::meas::algorithms::shapelet::Ellipse::doMeasure
bool doMeasure(const std::vector< PixelList > &pix, const std::vector< BVec > *psf, int order, int order2, int maxm, double sigma, long &flag, double thresh, DMatrix *cov=0)
Definition: Ellipse_meas.cc:39
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
bool lsst::meas::algorithms::shapelet::Ellipse::measureShapelet ( const std::vector< PixelList > &  pix,
BVec &  bret,
int  order,
int  order2,
int  maxm,
DMatrix *  bcov = 0 
) const

Definition at line 614 of file Ellipse.cc.

617  { return doMeasureShapelet(pix,0,b,order,order2,maxm,bCov); }
lsst::meas::algorithms::shapelet::Ellipse::doMeasureShapelet
bool doMeasureShapelet(const std::vector< PixelList > &pix, const std::vector< BVec > *psf, BVec &bret, int order, int order2, int maxm, DMatrix *bcov=0) const
Definition: Ellipse.cc:186
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
bool lsst::meas::algorithms::shapelet::Ellipse::measureShapelet ( const std::vector< PixelList > &  pix,
const std::vector< BVec > &  psf,
BVec &  bret,
int  order,
int  order2,
int  maxm,
DMatrix *  bcov = 0 
) const

Definition at line 608 of file Ellipse.cc.

612  { return doMeasureShapelet(pix,&psf,b,order,order2,maxm,bCov); }
lsst::meas::algorithms::shapelet::Ellipse::doMeasureShapelet
bool doMeasureShapelet(const std::vector< PixelList > &pix, const std::vector< BVec > *psf, BVec &bret, int order, int order2, int maxm, DMatrix *bcov=0) const
Definition: Ellipse.cc:186
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
Ellipse& lsst::meas::algorithms::shapelet::Ellipse::operator= ( const Ellipse &  e2)
inline

Definition at line 45 of file Ellipse.h.

46  {
47  _cen = e2.getCen();
48  _gamma = e2.getGamma();
49  _mu = std::complex<double>(e2.getMu(),e2.getTheta());
50  return *this;
51  }
void lsst::meas::algorithms::shapelet::Ellipse::peakCentroid ( const PixelList &  pix,
double  maxr 
)

Definition at line 354 of file CrudeMeasure.cc.

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  }
void lsst::meas::algorithms::shapelet::Ellipse::postShiftBy ( const std::complex< double > &  cen,
const std::complex< double > &  gamma,
const std::complex< double > &  mu 
)

Definition at line 135 of file Ellipse.cc.

139  {
140  // Current values are c1, g1, m1.
141  // The model is that a round galaxy is first transformed by (c1,g1,m1).
142  // And then by the provided values of (c2,g2,m2).
143  CalculateShift(_cen,_gamma,_mu,c2,g2,m2,_cen,_gamma,_mu);
144  }
lsst::meas::algorithms::shapelet::CalculateShift
void CalculateShift(std::complex< double > c1, std::complex< double > g1, std::complex< double > m1, std::complex< double > c2, std::complex< double > g2, std::complex< double > m2, std::complex< double > &c3, std::complex< double > &g3, std::complex< double > &m3)
Definition: Ellipse.cc:38
void lsst::meas::algorithms::shapelet::Ellipse::preShiftBy ( const std::complex< double > &  cen,
const std::complex< double > &  gamma,
const std::complex< double > &  mu 
)

Definition at line 124 of file Ellipse.cc.

128  {
129  // Current values are c2, g2, m2.
130  // The model is that a round galaxy is first transformed by (c1,g1,m1).
131  // And then by the current values of (c2,g2,m2).
132  CalculateShift(c1,g1,m1,_cen,_gamma,_mu,_cen,_gamma,_mu);
133  }
lsst::meas::algorithms::shapelet::CalculateShift
void CalculateShift(std::complex< double > c1, std::complex< double > g1, std::complex< double > m1, std::complex< double > c2, std::complex< double > g2, std::complex< double > m2, std::complex< double > &c3, std::complex< double > &g3, std::complex< double > &m3)
Definition: Ellipse.cc:38
void lsst::meas::algorithms::shapelet::Ellipse::removeRotation ( )

Definition at line 146 of file Ellipse.cc.

147  {
148  // This finds the rotation-free transformation that distorts a
149  // circle to the same final shape as the current transformation does.
150  // In otherwords, we move the rotation part from the end of the
151  // transformation to the beginning, since a rotation of a circle
152  // is a null operation.
153  //
154  // z' = exp(-m-it)/sqrt(1-|g|^2) ( z-c - g (z-c)* )
155  // = exp(-m)/sqrt(1-|g|^2) (exp(-it)(z-c) - g exp(-it) (z-c)*)
156  // = exp(-m)/sqrt(1-|g|^2) (exp(-it)(z-c) - g exp(-2it) (exp(-it)(z-c))*)
157  //
158  // So, g -> g exp(-2it)
159  // c -> c exp(-it)
160  // m -> m
161  // t -> 0
162 
163  if (imag(_mu) != 0.) {
164  std::complex<double> r = std::polar(1.,-imag(_mu));
165  _gamma *= r*r;
166  _cen *= r;
167  _mu = real(_mu);
168  }
169  }
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
void lsst::meas::algorithms::shapelet::Ellipse::setCen ( const std::complex< double > &  cen)
inline

Definition at line 137 of file Ellipse.h.

137 { _cen = cen; }
void lsst::meas::algorithms::shapelet::Ellipse::setGamma ( const std::complex< double > &  gamma)
inline

Definition at line 138 of file Ellipse.h.

138 { _gamma = gamma; }
void lsst::meas::algorithms::shapelet::Ellipse::setMu ( double  mu)
inline

Definition at line 139 of file Ellipse.h.

139 { _mu = std::complex<double>(mu,getTheta()); }
void lsst::meas::algorithms::shapelet::Ellipse::setTheta ( double  theta)
inline

Definition at line 140 of file Ellipse.h.

140 { _mu = std::complex<double>(getMu(),theta); }
void lsst::meas::algorithms::shapelet::Ellipse::unfixCen ( )
inline

Definition at line 146 of file Ellipse.h.

146 { _isFixedCen = false; }
void lsst::meas::algorithms::shapelet::Ellipse::unfixGam ( )
inline

Definition at line 147 of file Ellipse.h.

147 { _isFixedGamma = false; }
void lsst::meas::algorithms::shapelet::Ellipse::unfixMu ( )
inline

Definition at line 148 of file Ellipse.h.

148 { _isFixedMu = false; }
void lsst::meas::algorithms::shapelet::Ellipse::write ( std::ostream &  os) const
inline

Definition at line 156 of file Ellipse.h.

157  { os << _cen<<" "<<_gamma<<" "<<_mu; }

Member Data Documentation

std::complex<double> lsst::meas::algorithms::shapelet::Ellipse::_cen
private

Definition at line 176 of file Ellipse.h.

std::complex<double> lsst::meas::algorithms::shapelet::Ellipse::_gamma
private

Definition at line 177 of file Ellipse.h.

bool lsst::meas::algorithms::shapelet::Ellipse::_isFixedCen
private

Definition at line 180 of file Ellipse.h.

bool lsst::meas::algorithms::shapelet::Ellipse::_isFixedGamma
private

Definition at line 180 of file Ellipse.h.

bool lsst::meas::algorithms::shapelet::Ellipse::_isFixedMu
private

Definition at line 180 of file Ellipse.h.

std::complex<double> lsst::meas::algorithms::shapelet::Ellipse::_mu
private

Definition at line 178 of file Ellipse.h.

bool lsst::meas::algorithms::shapelet::Ellipse::_shouldDoTimings
private

Definition at line 182 of file Ellipse.h.

The documentation for this class was generated from the following files:
Generated on Thu Sep 24 2015 02:29:38 for LSSTApplications by   doxygen 1.8.5