LSSTApplications  8.0.0.0+107,8.0.0.1+13,9.1+18,9.2,master-g084aeec0a4,master-g0aced2eed8+6,master-g15627eb03c,master-g28afc54ef9,master-g3391ba5ea0,master-g3d0fb8ae5f,master-g4432ae2e89+36,master-g5c3c32f3ec+17,master-g60f1e072bb+1,master-g6a3ac32d1b,master-g76a88a4307+1,master-g7bce1f4e06+57,master-g8ff4092549+31,master-g98e65bf68e,master-ga6b77976b1+53,master-gae20e2b580+3,master-gb584cd3397+53,master-gc5448b162b+1,master-gc54cf9771d,master-gc69578ece6+1,master-gcbf758c456+22,master-gcec1da163f+63,master-gcf15f11bcc,master-gd167108223,master-gf44c96c709
LSSTDataManagementBasePackage
Classes | Functions | Variables
lsst::meas::algorithms::detail Namespace Reference

Classes

class  FftShifter
 
class  SdssShapeImpl
 

Functions

template<typename PixelT >
afwImage::Image< PixelT >::Ptr calcImageRealSpace (double const rad1, double const rad2, double const taperwidth)
 
std::pair< double, double > rotate (double x, double y, double angle)
 
template<typename PixelT >
afwImage::Image< PixelT >::Ptr calcImageKSpaceCplx (double const rad1, double const rad2, double const posAng, double const ellipticity)
 
template<typename PixelT >
afwImage::Image< PixelT >::Ptr calcImageKSpaceReal (double const rad1, double const rad2)
 
template<typename ImageT >
bool getAdaptiveMoments (ImageT const &mimage, double bkgd, double xcen, double ycen, double shiftmax, detail::SdssShapeImpl *shape, int maxIter, float tol1, float tol2)
 
template<typename ImageT >
std::pair< double, double > getFixedMomentsFlux (ImageT const &image, double bkgd, double xcen, double ycen, detail::SdssShapeImpl const &shape_)
 Return the flux of an object, using the aperture described by the SdssShape object. More...
 
template<typename PixelT >
lsst::afw::image::Image
< PixelT >::Ptr 		calcImageRealSpace (double const rad1, double const rad2, double const taper=0.1)
 
template<typename PixelT >
lsst::afw::image::Image
< PixelT >::Ptr 		calcImageKSpaceReal (double const rad1, double const rad2)
 
template<typename PixelT >
lsst::afw::image::Image
< PixelT >::Ptr 		calcImageKSpaceCplx (double const rad1, double const rad2, double const posAng, double const ell)
 

Variables

int const SDSS_SHAPE_MAX_ITER = 100
 
float const SDSS_SHAPE_TOL1 = 1.0e-5
 
float const SDSS_SHAPE_TOL2 = 1.0e-4
 

Function Documentation

template<typename PixelT >
lsst::afw::image::Image<PixelT>::Ptr lsst::meas::algorithms::detail::calcImageKSpaceCplx ( double const  rad1,
double const  rad2,
double const  posAng,
double const  ell 
)

todo

  • try sub pixel shift if it doesn't break even symmetry
  • put values directly in an Image
  • precompute the plan

Definition at line 445 of file SincFlux.cc.

447  {
448 
449  // we only need a half-width due to symmetry
450  // make the hwid 2*rad2 so we have some buffer space and round up to the next power of 2
451  int log2 = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
452  if (log2 < 3) { log2 = 3; }
453  int hwid = pow(2, log2);
454  int wid = 2*hwid - 1;
455  int xcen = wid/2, ycen = wid/2;
456  FftShifter fftshift(wid);
457 
458  boost::shared_array<std::complex<double> > cimg(new std::complex<double>[wid*wid]);
459  std::complex<double> *c = cimg.get();
460  // fftplan args: nx, ny, *in, *out, direction, flags
461  // - done in-situ if *in == *out
462  fftw_plan plan = fftw_plan_dft_2d(wid, wid,
463  reinterpret_cast<fftw_complex*>(c),
464  reinterpret_cast<fftw_complex*>(c),
465  FFTW_BACKWARD, FFTW_ESTIMATE);
466 
467  // compute the k-space values and put them in the cimg array
468  double const twoPiRad1 = afwGeom::TWOPI*rad1;
469  double const twoPiRad2 = afwGeom::TWOPI*rad2;
470  double const scale = (1.0 - ellipticity);
471  for (int iY = 0; iY < wid; ++iY) {
472  int const fY = fftshift.shift(iY);
473  double const ky = (static_cast<double>(iY) - ycen)/wid;
474 
475  for (int iX = 0; iX < wid; ++iX) {
476 
477  int const fX = fftshift.shift(iX);
478  double const kx = static_cast<double>(iX - xcen)/wid;
479 
480  // rotate
481  std::pair<double, double> coo = rotate(kx, ky, posAng);
482  double kxr = coo.first;
483  double kyr = coo.second;
484  // rescale
485  double const k = ::sqrt(kxr*kxr + scale*scale*kyr*kyr);
486 
487  double const airy1 = (rad1 > 0 ? rad1*J1(twoPiRad1*k) : 0.0)/k;
488  double const airy2 = rad2*J1(twoPiRad2*k)/k;
489  double const airy = airy2 - airy1;
490 
491  c[fY*wid + fX] = std::complex<double>(scale*airy, 0.0);
492  }
493  }
494  c[0] = scale*afwGeom::PI*(rad2*rad2 - rad1*rad1);
495 
496  // perform the fft and clean up after ourselves
497  fftw_execute(plan);
498  fftw_destroy_plan(plan);
499 
500  // put the coefficients into an image
501  typename afwImage::Image<PixelT>::Ptr coeffImage =
502  boost::make_shared<afwImage::Image<PixelT> >(afwGeom::ExtentI(wid, wid), 0.0);
503 
504  for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
505  int iX = 0;
506  typename afwImage::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
507  for (typename afwImage::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY); ptr != end; ++ptr) {
508  int fX = fftshift.shift(iX);
509  int fY = fftshift.shift(iY);
510  *ptr = static_cast<PixelT>(c[fY*wid + fX].real()/(wid*wid));
511  iX++;
512  }
513  }
514 
515  // reset the origin to be the middle of the image
516  coeffImage->setXY0(-wid/2, -wid/2);
517  return coeffImage;
518  }
lsst::afw::geom::PI
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
lsst::afw::image::ImageBase::row_end
x_iterator row_end(int y) const
Return an x_iterator to the end of the y&#39;th row.
Definition: Image.h:324
lsst::meas::algorithms::detail::rotate
std::pair< double, double > rotate(double x, double y, double angle)
Definition: SincFlux.cc:432
lsst::afw::geom::ExtentI
Extent< int, 2 > ExtentI
Definition: Extent.h:347
lsst::afw::geom::TWOPI
double const TWOPI
Definition: Angle.h:19
lsst::afw::image::ImageBase::row_begin
x_iterator row_begin(int y) const
Definition: Image.h:319
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:415
template<typename PixelT >
template lsst::afw::image::Image< double >::Ptr lsst::meas::algorithms::detail::calcImageKSpaceCplx< double > ( double const  rad1,
double const  rad2,
double const  posAng,
double const  ellipticity 
)

todo

  • try sub pixel shift if it doesn't break even symmetry
  • put values directly in an Image
  • precompute the plan

Definition at line 445 of file SincFlux.cc.

447  {
448 
449  // we only need a half-width due to symmetry
450  // make the hwid 2*rad2 so we have some buffer space and round up to the next power of 2
451  int log2 = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
452  if (log2 < 3) { log2 = 3; }
453  int hwid = pow(2, log2);
454  int wid = 2*hwid - 1;
455  int xcen = wid/2, ycen = wid/2;
456  FftShifter fftshift(wid);
457 
458  boost::shared_array<std::complex<double> > cimg(new std::complex<double>[wid*wid]);
459  std::complex<double> *c = cimg.get();
460  // fftplan args: nx, ny, *in, *out, direction, flags
461  // - done in-situ if *in == *out
462  fftw_plan plan = fftw_plan_dft_2d(wid, wid,
463  reinterpret_cast<fftw_complex*>(c),
464  reinterpret_cast<fftw_complex*>(c),
465  FFTW_BACKWARD, FFTW_ESTIMATE);
466 
467  // compute the k-space values and put them in the cimg array
468  double const twoPiRad1 = afwGeom::TWOPI*rad1;
469  double const twoPiRad2 = afwGeom::TWOPI*rad2;
470  double const scale = (1.0 - ellipticity);
471  for (int iY = 0; iY < wid; ++iY) {
472  int const fY = fftshift.shift(iY);
473  double const ky = (static_cast<double>(iY) - ycen)/wid;
474 
475  for (int iX = 0; iX < wid; ++iX) {
476 
477  int const fX = fftshift.shift(iX);
478  double const kx = static_cast<double>(iX - xcen)/wid;
479 
480  // rotate
481  std::pair<double, double> coo = rotate(kx, ky, posAng);
482  double kxr = coo.first;
483  double kyr = coo.second;
484  // rescale
485  double const k = ::sqrt(kxr*kxr + scale*scale*kyr*kyr);
486 
487  double const airy1 = (rad1 > 0 ? rad1*J1(twoPiRad1*k) : 0.0)/k;
488  double const airy2 = rad2*J1(twoPiRad2*k)/k;
489  double const airy = airy2 - airy1;
490 
491  c[fY*wid + fX] = std::complex<double>(scale*airy, 0.0);
492  }
493  }
494  c[0] = scale*afwGeom::PI*(rad2*rad2 - rad1*rad1);
495 
496  // perform the fft and clean up after ourselves
497  fftw_execute(plan);
498  fftw_destroy_plan(plan);
499 
500  // put the coefficients into an image
501  typename afwImage::Image<PixelT>::Ptr coeffImage =
502  boost::make_shared<afwImage::Image<PixelT> >(afwGeom::ExtentI(wid, wid), 0.0);
503 
504  for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
505  int iX = 0;
506  typename afwImage::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
507  for (typename afwImage::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY); ptr != end; ++ptr) {
508  int fX = fftshift.shift(iX);
509  int fY = fftshift.shift(iY);
510  *ptr = static_cast<PixelT>(c[fY*wid + fX].real()/(wid*wid));
511  iX++;
512  }
513  }
514 
515  // reset the origin to be the middle of the image
516  coeffImage->setXY0(-wid/2, -wid/2);
517  return coeffImage;
518  }
lsst::afw::geom::PI
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
lsst::afw::image::ImageBase::row_end
x_iterator row_end(int y) const
Return an x_iterator to the end of the y&#39;th row.
Definition: Image.h:324
lsst::meas::algorithms::detail::rotate
std::pair< double, double > rotate(double x, double y, double angle)
Definition: SincFlux.cc:432
lsst::afw::geom::ExtentI
Extent< int, 2 > ExtentI
Definition: Extent.h:347
lsst::afw::geom::TWOPI
double const TWOPI
Definition: Angle.h:19
lsst::afw::image::ImageBase::row_begin
x_iterator row_begin(int y) const
Definition: Image.h:319
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:415
template<typename PixelT >
lsst::afw::image::Image<PixelT>::Ptr lsst::meas::algorithms::detail::calcImageKSpaceReal ( double const  rad1,
double const  rad2 
)

Definition at line 527 of file SincFlux.cc.

527  {
528 
529  // we only need a half-width due to symmertry
530  // make the hwid 2*rad2 so we have some buffer space and round up to the next power of 2
531  int log2 = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
532  if (log2 < 3) { log2 = 3; }
533  int hwid = pow(2, log2);
534  int wid = 2*hwid - 1;
535  int xcen = wid/2, ycen = wid/2;
536  FftShifter fftshift(wid);
537 
538  boost::shared_array<double> cimg(new double[wid*wid]);
539  double *c = cimg.get();
540  // fftplan args: nx, ny, *in, *out, kindx, kindy, flags
541  // - done in-situ if *in == *out
542  fftw_plan plan = fftw_plan_r2r_2d(wid, wid, c, c, FFTW_R2HC, FFTW_R2HC, FFTW_ESTIMATE);
543 
544  // compute the k-space values and put them in the cimg array
545  double const twoPiRad1 = afwGeom::TWOPI*rad1;
546  double const twoPiRad2 = afwGeom::TWOPI*rad2;
547  for (int iY = 0; iY < wid; ++iY) {
548 
549  int const fY = fftshift.shift(iY);
550  double const ky = (static_cast<double>(iY) - ycen)/wid;
551 
552  for (int iX = 0; iX < wid; ++iX) {
553  int const fX = fftshift.shift(iX);
554 
555  // emacs indent breaks if this isn't separte
556  double const iXcen = static_cast<double>(iX - xcen);
557  double const kx = iXcen/wid;
558 
559  double const k = ::sqrt(kx*kx + ky*ky);
560  double const airy1 = (rad1 > 0 ? rad1*J1(twoPiRad1*k) : 0.0)/k;
561  double const airy2 = rad2*J1(twoPiRad2*k)/k;
562  double const airy = airy2 - airy1;
563  c[fY*wid + fX] = airy;
564 
565  }
566  }
567  int fxy = fftshift.shift(wid/2);
568  c[fxy*wid + fxy] = afwGeom::PI*(rad2*rad2 - rad1*rad1);
569 
570  // perform the fft and clean up after ourselves
571  fftw_execute(plan);
572  fftw_destroy_plan(plan);
573 
574  // put the coefficients into an image
575  typename afwImage::Image<PixelT>::Ptr coeffImage =
576  boost::make_shared<afwImage::Image<PixelT> >(afwGeom::ExtentI(wid, wid), 0.0);
577 
578  for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
579  int iX = 0;
580  typename afwImage::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
581  for (typename afwImage::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY); ptr != end; ++ptr) {
582 
583  // now need to reflect the quadrant we solved to the other three
584  int fX = iX < hwid ? hwid - iX - 1 : iX - hwid + 1;
585  int fY = iY < hwid ? hwid - iY - 1 : iY - hwid + 1;
586  *ptr = static_cast<PixelT>(c[fY*wid + fX]/(wid*wid));
587  iX++;
588  }
589  }
590 
591  // reset the origin to be the middle of the image
592  coeffImage->setXY0(-wid/2, -wid/2);
593  return coeffImage;
594  }
lsst::afw::geom::PI
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
lsst::afw::image::ImageBase::row_end
x_iterator row_end(int y) const
Return an x_iterator to the end of the y&#39;th row.
Definition: Image.h:324
lsst::afw::geom::ExtentI
Extent< int, 2 > ExtentI
Definition: Extent.h:347
lsst::afw::geom::TWOPI
double const TWOPI
Definition: Angle.h:19
lsst::afw::image::ImageBase::row_begin
x_iterator row_begin(int y) const
Definition: Image.h:319
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:415
template<typename PixelT >
template lsst::afw::image::Image< double >::Ptr lsst::meas::algorithms::detail::calcImageKSpaceReal< double > ( double const  rad1,
double const  rad2 
)

Definition at line 527 of file SincFlux.cc.

527  {
528 
529  // we only need a half-width due to symmertry
530  // make the hwid 2*rad2 so we have some buffer space and round up to the next power of 2
531  int log2 = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
532  if (log2 < 3) { log2 = 3; }
533  int hwid = pow(2, log2);
534  int wid = 2*hwid - 1;
535  int xcen = wid/2, ycen = wid/2;
536  FftShifter fftshift(wid);
537 
538  boost::shared_array<double> cimg(new double[wid*wid]);
539  double *c = cimg.get();
540  // fftplan args: nx, ny, *in, *out, kindx, kindy, flags
541  // - done in-situ if *in == *out
542  fftw_plan plan = fftw_plan_r2r_2d(wid, wid, c, c, FFTW_R2HC, FFTW_R2HC, FFTW_ESTIMATE);
543 
544  // compute the k-space values and put them in the cimg array
545  double const twoPiRad1 = afwGeom::TWOPI*rad1;
546  double const twoPiRad2 = afwGeom::TWOPI*rad2;
547  for (int iY = 0; iY < wid; ++iY) {
548 
549  int const fY = fftshift.shift(iY);
550  double const ky = (static_cast<double>(iY) - ycen)/wid;
551 
552  for (int iX = 0; iX < wid; ++iX) {
553  int const fX = fftshift.shift(iX);
554 
555  // emacs indent breaks if this isn't separte
556  double const iXcen = static_cast<double>(iX - xcen);
557  double const kx = iXcen/wid;
558 
559  double const k = ::sqrt(kx*kx + ky*ky);
560  double const airy1 = (rad1 > 0 ? rad1*J1(twoPiRad1*k) : 0.0)/k;
561  double const airy2 = rad2*J1(twoPiRad2*k)/k;
562  double const airy = airy2 - airy1;
563  c[fY*wid + fX] = airy;
564 
565  }
566  }
567  int fxy = fftshift.shift(wid/2);
568  c[fxy*wid + fxy] = afwGeom::PI*(rad2*rad2 - rad1*rad1);
569 
570  // perform the fft and clean up after ourselves
571  fftw_execute(plan);
572  fftw_destroy_plan(plan);
573 
574  // put the coefficients into an image
575  typename afwImage::Image<PixelT>::Ptr coeffImage =
576  boost::make_shared<afwImage::Image<PixelT> >(afwGeom::ExtentI(wid, wid), 0.0);
577 
578  for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
579  int iX = 0;
580  typename afwImage::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
581  for (typename afwImage::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY); ptr != end; ++ptr) {
582 
583  // now need to reflect the quadrant we solved to the other three
584  int fX = iX < hwid ? hwid - iX - 1 : iX - hwid + 1;
585  int fY = iY < hwid ? hwid - iY - 1 : iY - hwid + 1;
586  *ptr = static_cast<PixelT>(c[fY*wid + fX]/(wid*wid));
587  iX++;
588  }
589  }
590 
591  // reset the origin to be the middle of the image
592  coeffImage->setXY0(-wid/2, -wid/2);
593  return coeffImage;
594  }
lsst::afw::geom::PI
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
lsst::afw::image::ImageBase::row_end
x_iterator row_end(int y) const
Return an x_iterator to the end of the y&#39;th row.
Definition: Image.h:324
lsst::afw::geom::ExtentI
Extent< int, 2 > ExtentI
Definition: Extent.h:347
lsst::afw::geom::TWOPI
double const TWOPI
Definition: Angle.h:19
lsst::afw::image::ImageBase::row_begin
x_iterator row_begin(int y) const
Definition: Image.h:319
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:415
template<typename PixelT >
lsst::afw::image::Image<PixelT>::Ptr lsst::meas::algorithms::detail::calcImageRealSpace ( double const  rad1,
double const  rad2,
double const  taper = 0.1 
)

Definition at line 326 of file SincFlux.cc.

326  {
327 
328  PixelT initweight = 0.0; // initialize the coeff values
329 
330  int log2 = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
331  if (log2 < 3) { log2 = 3; }
332  int hwid = pow(2, log2);
333  int width = 2*hwid - 1;
334 
335  int const xwidth = width;
336  int const ywidth = width;
337 
338  int const x0 = -xwidth/2;
339  int const y0 = -ywidth/2;
340 
341  // create an image to hold the coefficient image
342  typename afwImage::Image<PixelT>::Ptr coeffImage =
343  boost::make_shared<afwImage::Image<PixelT> >(afwGeom::ExtentI(xwidth, ywidth), initweight);
344  coeffImage->setXY0(x0, y0);
345 
346  // create the aperture function object
347  // determine the radius to use that makes 'radius' the effective radius of the aperture
348  double tolerance = 1.0e-12;
349  double dr = 1.0e-6;
350  double err = 2.0*tolerance;
351  double apEff = afwGeom::PI*rad2*rad2;
352  double radIn = rad2;
353  int maxIt = 20;
354  int i = 0;
355  while (err > tolerance && i < maxIt) {
356  CircApPolar<double> apPolar1(radIn, taperwidth);
357  CircApPolar<double> apPolar2(radIn+dr, taperwidth);
358  double a1 = afwGeom::TWOPI * afwMath::integrate(apPolar1, 0.0, radIn+taperwidth, tolerance);
359  double a2 = afwGeom::TWOPI * afwMath::integrate(apPolar2, 0.0, radIn+dr+taperwidth, tolerance);
360  double dadr = (a2 - a1)/dr;
361  double radNew = radIn - (a1 - apEff)/dadr;
362  err = (a1 - apEff)/apEff;
363  radIn = radNew;
364  i++;
365  }
366  CircularAperture<double> ap(rad1, rad2, taperwidth);
367 
368 
369  /* ******************************* */
370  // integrate over the aperture
371 
372  // the limits of the integration over the sinc aperture
373  double const limit = rad2 + taperwidth;
374  double const x1 = -limit;
375  double const x2 = limit;
376  double const y1 = -limit;
377  double const y2 = limit;
378 
379  for (int iY = y0; iY != y0 + coeffImage->getHeight(); ++iY) {
380  int iX = x0;
381  typename afwImage::Image<PixelT>::x_iterator end = coeffImage->row_end(iY-y0);
382  for (typename afwImage::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY-y0); ptr != end; ++ptr) {
383 
384  // create a sinc function in the CircularAperture at our location
385  SincAperture<double> sincAp(ap, iX, iY);
386 
387  // integrate the sinc
388  PixelT integral = afwMath::integrate2d(sincAp, x1, x2, y1, y2, 1.0e-8);
389 
390  // we actually integrated function+1.0 and now must subtract the excess volume
391  // - just force it to zero in the corners
392  double const dx = iX;
393  double const dy = iY;
394  *ptr = (std::sqrt(dx*dx + dy*dy) < xwidth/2) ?
395  integral - (x2 - x1)*(y2 - y1) : 0.0;
396  ++iX;
397  }
398  }
399 
400 
401  double sum = 0.0;
402  for (int iY = y0; iY != y0 + coeffImage->getHeight(); ++iY) {
403  typename afwImage::Image<PixelT>::x_iterator end = coeffImage->row_end(iY-y0);
404  for (typename afwImage::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY-y0); ptr != end; ++ptr) {
405  sum += *ptr;
406  }
407  }
408 
409 #if 0 // debugging
410  coeffImage->writeFits("cimage.fits");
411 #endif
412 
413  return coeffImage;
414  }
lsst::afw::image::ImageBase::setXY0
void setXY0(geom::Point2I const origin)
Definition: Image.h:362
dx
double dx
Definition: ImageUtils.cc:90
lsst::afw::math::integrate2d
BinaryFunctionT::result_type integrate2d(BinaryFunctionT func, typename BinaryFunctionT::first_argument_type const x1, typename BinaryFunctionT::first_argument_type const x2, typename BinaryFunctionT::second_argument_type const y1, typename BinaryFunctionT::second_argument_type const y2, double eps=1.0e-6)
The 2D integrator.
Definition: Integrate.h:977
dy
double dy
Definition: ImageUtils.cc:90
lsst::afw::image::Image::writeFits
void writeFits(std::string const &fileName, boost::shared_ptr< lsst::daf::base::PropertySet const > metadata=boost::shared_ptr< lsst::daf::base::PropertySet const >(), std::string const &mode="w") const
Write an image to a regular FITS file.
lsst::afw::geom::PI
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
lsst::afw::image::ImageBase::row_end
x_iterator row_end(int y) const
Return an x_iterator to the end of the y&#39;th row.
Definition: Image.h:324
lsst::afw::image::ImageBase::getHeight
int getHeight() const
Return the number of rows in the image.
Definition: Image.h:239
sum
double sum
Definition: NaiveFlux.cc:137
lsst::afw::geom::ExtentI
Extent< int, 2 > ExtentI
Definition: Extent.h:347
lsst::afw::geom::TWOPI
double const TWOPI
Definition: Angle.h:19
lsst::afw::image::ImageBase::row_begin
x_iterator row_begin(int y) const
Definition: Image.h:319
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:415
lsst::afw::math::integrate
UnaryFunctionT::result_type integrate(UnaryFunctionT func, typename UnaryFunctionT::argument_type const a, typename UnaryFunctionT::argument_type const b, double eps=1.0e-6)
The 1D integrator.
Definition: Integrate.h:917
template<typename PixelT >
template lsst::afw::image::Image< double >::Ptr lsst::meas::algorithms::detail::calcImageRealSpace< double > ( double const  rad1,
double const  rad2,
double const  taperwidth 
)

Definition at line 326 of file SincFlux.cc.

326  {
327 
328  PixelT initweight = 0.0; // initialize the coeff values
329 
330  int log2 = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
331  if (log2 < 3) { log2 = 3; }
332  int hwid = pow(2, log2);
333  int width = 2*hwid - 1;
334 
335  int const xwidth = width;
336  int const ywidth = width;
337 
338  int const x0 = -xwidth/2;
339  int const y0 = -ywidth/2;
340 
341  // create an image to hold the coefficient image
342  typename afwImage::Image<PixelT>::Ptr coeffImage =
343  boost::make_shared<afwImage::Image<PixelT> >(afwGeom::ExtentI(xwidth, ywidth), initweight);
344  coeffImage->setXY0(x0, y0);
345 
346  // create the aperture function object
347  // determine the radius to use that makes 'radius' the effective radius of the aperture
348  double tolerance = 1.0e-12;
349  double dr = 1.0e-6;
350  double err = 2.0*tolerance;
351  double apEff = afwGeom::PI*rad2*rad2;
352  double radIn = rad2;
353  int maxIt = 20;
354  int i = 0;
355  while (err > tolerance && i < maxIt) {
356  CircApPolar<double> apPolar1(radIn, taperwidth);
357  CircApPolar<double> apPolar2(radIn+dr, taperwidth);
358  double a1 = afwGeom::TWOPI * afwMath::integrate(apPolar1, 0.0, radIn+taperwidth, tolerance);
359  double a2 = afwGeom::TWOPI * afwMath::integrate(apPolar2, 0.0, radIn+dr+taperwidth, tolerance);
360  double dadr = (a2 - a1)/dr;
361  double radNew = radIn - (a1 - apEff)/dadr;
362  err = (a1 - apEff)/apEff;
363  radIn = radNew;
364  i++;
365  }
366  CircularAperture<double> ap(rad1, rad2, taperwidth);
367 
368 
369  /* ******************************* */
370  // integrate over the aperture
371 
372  // the limits of the integration over the sinc aperture
373  double const limit = rad2 + taperwidth;
374  double const x1 = -limit;
375  double const x2 = limit;
376  double const y1 = -limit;
377  double const y2 = limit;
378 
379  for (int iY = y0; iY != y0 + coeffImage->getHeight(); ++iY) {
380  int iX = x0;
381  typename afwImage::Image<PixelT>::x_iterator end = coeffImage->row_end(iY-y0);
382  for (typename afwImage::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY-y0); ptr != end; ++ptr) {
383 
384  // create a sinc function in the CircularAperture at our location
385  SincAperture<double> sincAp(ap, iX, iY);
386 
387  // integrate the sinc
388  PixelT integral = afwMath::integrate2d(sincAp, x1, x2, y1, y2, 1.0e-8);
389 
390  // we actually integrated function+1.0 and now must subtract the excess volume
391  // - just force it to zero in the corners
392  double const dx = iX;
393  double const dy = iY;
394  *ptr = (std::sqrt(dx*dx + dy*dy) < xwidth/2) ?
395  integral - (x2 - x1)*(y2 - y1) : 0.0;
396  ++iX;
397  }
398  }
399 
400 
401  double sum = 0.0;
402  for (int iY = y0; iY != y0 + coeffImage->getHeight(); ++iY) {
403  typename afwImage::Image<PixelT>::x_iterator end = coeffImage->row_end(iY-y0);
404  for (typename afwImage::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY-y0); ptr != end; ++ptr) {
405  sum += *ptr;
406  }
407  }
408 
409 #if 0 // debugging
410  coeffImage->writeFits("cimage.fits");
411 #endif
412 
413  return coeffImage;
414  }
lsst::afw::image::ImageBase::setXY0
void setXY0(geom::Point2I const origin)
Definition: Image.h:362
dx
double dx
Definition: ImageUtils.cc:90
lsst::afw::math::integrate2d
BinaryFunctionT::result_type integrate2d(BinaryFunctionT func, typename BinaryFunctionT::first_argument_type const x1, typename BinaryFunctionT::first_argument_type const x2, typename BinaryFunctionT::second_argument_type const y1, typename BinaryFunctionT::second_argument_type const y2, double eps=1.0e-6)
The 2D integrator.
Definition: Integrate.h:977
dy
double dy
Definition: ImageUtils.cc:90
lsst::afw::image::Image::writeFits
void writeFits(std::string const &fileName, boost::shared_ptr< lsst::daf::base::PropertySet const > metadata=boost::shared_ptr< lsst::daf::base::PropertySet const >(), std::string const &mode="w") const
Write an image to a regular FITS file.
lsst::afw::geom::PI
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
lsst::afw::image::ImageBase::row_end
x_iterator row_end(int y) const
Return an x_iterator to the end of the y&#39;th row.
Definition: Image.h:324
lsst::afw::image::ImageBase::getHeight
int getHeight() const
Return the number of rows in the image.
Definition: Image.h:239
sum
double sum
Definition: NaiveFlux.cc:137
lsst::afw::geom::ExtentI
Extent< int, 2 > ExtentI
Definition: Extent.h:347
lsst::afw::geom::TWOPI
double const TWOPI
Definition: Angle.h:19
lsst::afw::image::ImageBase::row_begin
x_iterator row_begin(int y) const
Definition: Image.h:319
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:415
lsst::afw::math::integrate
UnaryFunctionT::result_type integrate(UnaryFunctionT func, typename UnaryFunctionT::argument_type const a, typename UnaryFunctionT::argument_type const b, double eps=1.0e-6)
The 1D integrator.
Definition: Integrate.h:917
template<typename ImageT >
bool lsst::meas::algorithms::detail::getAdaptiveMoments ( ImageT const &  mimage,
double  bkgd,
double  xcen,
double  ycen,
double  shiftmax,
detail::SdssShapeImpl *  shape,
int  maxIter,
float  tol1,
float  tol2 
)

Workhorse for adaptive moments

Calculate adaptive moments from an image

The moments are measured iteratively with a Gaussian window with width equal to the second moments from the previous iteration.

Parameters
mimagethe data to process
bkgdbackground level
xcenx-centre of object
yceny-centre of object
shiftmaxmax allowed centroid shift
shapea place to store desired data
maxIterMaximum number of iterations
tol1Convergence tolerance for e1,e2
tol2Convergence tolerance for FWHM

Definition at line 453 of file SdssShape.cc.

455 {
456  double I0 = 0; // amplitude of best-fit Gaussian
457  double sum; // sum of intensity*weight
458  double sumx, sumy; // sum ((int)[xy])*intensity*weight
459  double sumxx, sumxy, sumyy; // sum {x^2,xy,y^2}*intensity*weight
460  double sums4; // sum intensity*weight*exponent^2
461  float const xcen0 = xcen; // initial centre
462  float const ycen0 = ycen; // of object
463 
464  double sigma11W = 1.5; // quadratic moments of the
465  double sigma12W = 0.0; // weighting fcn;
466  double sigma22W = 1.5; // xx, xy, and yy
467 
468  double w11 = -1, w12 = -1, w22 = -1; // current weights for moments; always set when iter == 0
469  float e1_old = 1e6, e2_old = 1e6; // old values of shape parameters e1 and e2
470  float sigma11_ow_old = 1e6; // previous version of sigma11_ow
471 
472  typename ImageAdaptor<ImageT>::Image const &image = ImageAdaptor<ImageT>().getImage(mimage);
473 
474  if (lsst::utils::isnan(xcen) || lsst::utils::isnan(ycen)) {
475  // Can't do anything
476  shape->setFlag(SdssShapeImpl::UNWEIGHTED_BAD);
477  return false;
478  }
479 
480  bool interpflag = false; // interpolate finer than a pixel?
481  lsst::afw::geom::BoxI bbox;
482  int iter = 0; // iteration number
483  for (; iter < maxIter; iter++) {
484  bbox = set_amom_bbox(image.getWidth(), image.getHeight(), xcen, ycen, sigma11W, sigma12W, sigma22W);
485 
486  boost::tuple<std::pair<bool, double>, double, double, double> weights =
487  getWeights(sigma11W, sigma12W, sigma22W);
488  if (!weights.get<0>().first) {
489  shape->setFlag(SdssShapeImpl::UNWEIGHTED);
490  break;
491  }
492 
493  double const detW = weights.get<0>().second;
494 
495 #if 0 // this form was numerically unstable on my G4 powerbook
496  assert(detW >= 0.0);
497 #else
498  assert(sigma11W*sigma22W >= sigma12W*sigma12W - std::numeric_limits<float>::epsilon());
499 #endif
500 
501  {
502  const double ow11 = w11; // old
503  const double ow12 = w12; // values
504  const double ow22 = w22; // of w11, w12, w22
505 
506  w11 = weights.get<1>();
507  w12 = weights.get<2>();
508  w22 = weights.get<3>();
509 
510  if (shouldInterp(sigma11W, sigma22W, detW)) {
511  if (!interpflag) {
512  interpflag = true; // N.b.: stays set for this object
513  if (iter > 0) {
514  sigma11_ow_old = 1.e6; // force at least one more iteration
515  w11 = ow11;
516  w12 = ow12;
517  w22 = ow22;
518  iter--; // we didn't update wXX
519  }
520  }
521  }
522  }
523 
524  if (calcmom<false>(image, xcen, ycen, bbox, bkgd, interpflag, w11, w12, w22,
525  &I0, &sum, &sumx, &sumy, &sumxx, &sumxy, &sumyy, &sums4) < 0) {
526  shape->setFlag(SdssShapeImpl::UNWEIGHTED);
527  break;
528  }
529 #if 0
530 /*
531  * Find new centre
532  *
533  * This is only needed if we update the centre; if we use the input position we've already done the work
534  */
535  xcen = sumx/sum;
536  ycen = sumy/sum;
537 #endif
538  shape->setX(sumx/sum); // update centroid. N.b. we're not setting errors here
539  shape->setY(sumy/sum);
540 
541  if (fabs(shape->getX() - xcen0) > shiftmax || fabs(shape->getY() - ycen0) > shiftmax) {
542  shape->setFlag(SdssShapeImpl::SHIFT);
543  }
544 /*
545  * OK, we have the centre. Proceed to find the second moments.
546  */
547  float const sigma11_ow = sumxx/sum; // quadratic moments of
548  float const sigma22_ow = sumyy/sum; // weight*object
549  float const sigma12_ow = sumxy/sum; // xx, xy, and yy
550 
551  if (sigma11_ow <= 0 || sigma22_ow <= 0) {
552  shape->setFlag(SdssShapeImpl::UNWEIGHTED);
553  break;
554  }
555 
556  float const d = sigma11_ow + sigma22_ow; // current values of shape parameters
557  float const e1 = (sigma11_ow - sigma22_ow)/d;
558  float const e2 = 2.0*sigma12_ow/d;
559 /*
560  * Did we converge?
561  */
562  if (iter > 0 &&
563  fabs(e1 - e1_old) < tol1 && fabs(e2 - e2_old) < tol1 &&
564  fabs(sigma11_ow/sigma11_ow_old - 1.0) < tol2 ) {
565  break; // yes; we converged
566  }
567 
568  e1_old = e1;
569  e2_old = e2;
570  sigma11_ow_old = sigma11_ow;
571 /*
572  * Didn't converge, calculate new values for weighting function
573  *
574  * The product of two Gaussians is a Gaussian:
575  * <x^2 exp(-a x^2 - 2bxy - cy^2) exp(-Ax^2 - 2Bxy - Cy^2)> =
576  * <x^2 exp(-(a + A) x^2 - 2(b + B)xy - (c + C)y^2)>
577  * i.e. the inverses of the covariances matrices add.
578  *
579  * We know sigmaXX_ow and sigmaXXW, the covariances of the weighted object
580  * and of the weights themselves. We can estimate the object's covariance as
581  * sigmaXX_ow^-1 - sigmaXXW^-1
582  * and, as we want to find a set of weights with the _same_ covariance as the
583  * object we take this to be the an estimate of our correct weights.
584  *
585  * N.b. This assumes that the object is roughly Gaussian.
586  * Consider the object:
587  * O == delta(x + p) + delta(x - p)
588  * the covariance of the weighted object is equal to that of the unweighted
589  * object, and this prescription fails badly. If we detect this, we set
590  * the UNWEIGHTED flag, and calculate the UNweighted moments
591  * instead.
592  */
593  {
594  float n11, n12, n22; // elements of inverse of next guess at weighting function
595  float ow11, ow12, ow22; // elements of inverse of sigmaXX_ow
596 
597  boost::tuple<std::pair<bool, double>, double, double, double> weights =
598  getWeights(sigma11_ow, sigma12_ow, sigma22_ow);
599  if (!weights.get<0>().first) {
600  shape->setFlag(SdssShapeImpl::UNWEIGHTED);
601  break;
602  }
603 
604  ow11 = weights.get<1>();
605  ow12 = weights.get<2>();
606  ow22 = weights.get<3>();
607 
608  n11 = ow11 - w11;
609  n12 = ow12 - w12;
610  n22 = ow22 - w22;
611 
612  weights = getWeights(n11, n12, n22, false);
613  if (!weights.get<0>().first) {
614  // product-of-Gaussians assumption failed
615  shape->setFlag(SdssShapeImpl::UNWEIGHTED);
616  break;
617  }
618 
619  sigma11W = weights.get<1>();
620  sigma12W = weights.get<2>();
621  sigma22W = weights.get<3>();
622  }
623 
624  if (sigma11W <= 0 || sigma22W <= 0) {
625  shape->setFlag(SdssShapeImpl::UNWEIGHTED);
626  break;
627  }
628  }
629 
630  if (iter == maxIter) {
631  shape->setFlag(SdssShapeImpl::UNWEIGHTED);
632  shape->setFlag(SdssShapeImpl::MAXITER);
633  }
634 
635  if (sumxx + sumyy == 0.0) {
636  shape->setFlag(SdssShapeImpl::UNWEIGHTED);
637  }
638 /*
639  * Problems; try calculating the un-weighted moments
640  */
641  if (shape->getFlag(SdssShapeImpl::UNWEIGHTED)) {
642  w11 = w22 = w12 = 0;
643  if (calcmom<false>(image, xcen, ycen, bbox, bkgd, interpflag, w11, w12, w22,
644  &I0, &sum, &sumx, &sumy, &sumxx, &sumxy, &sumyy, NULL) < 0 || sum <= 0) {
645  shape->resetFlag(SdssShapeImpl::UNWEIGHTED);
646  shape->setFlag(SdssShapeImpl::UNWEIGHTED_BAD);
647 
648  if (sum > 0) {
649  shape->setIxx(1/12.0); // a single pixel
650  shape->setIxy(0.0);
651  shape->setIyy(1/12.0);
652  }
653 
654  return false;
655  }
656 
657  sigma11W = sumxx/sum; // estimate of object moments
658  sigma12W = sumxy/sum; // usually, object == weight
659  sigma22W = sumyy/sum; // at this point
660  }
661 
662  shape->setI0(I0);
663  shape->setIxx(sigma11W);
664  shape->setIxy(sigma12W);
665  shape->setIyy(sigma22W);
666  shape->setIxy4(sums4/sum);
667 
668  if (shape->getIxx() + shape->getIyy() != 0.0) {
669  int const ix = lsst::afw::image::positionToIndex(xcen);
670  int const iy = lsst::afw::image::positionToIndex(ycen);
671 
672  if (ix >= 0 && ix < mimage.getWidth() && iy >= 0 && iy < mimage.getHeight()) {
673  float const bkgd_var =
674  ImageAdaptor<ImageT>().getVariance(mimage, ix, iy); // XXX Overestimate as it includes object
675 
676  if (bkgd_var > 0.0) { // NaN is not > 0.0
677  if (!(shape->getFlag(SdssShapeImpl::UNWEIGHTED))) {
678  detail::SdssShapeImpl::Matrix4 fisher = calc_fisher(*shape, bkgd_var); // Fisher matrix
679  shape->setCovar(fisher.inverse());
680  }
681  }
682  }
683  }
684 
685  return true;
686 }
lsst::afw::image::positionToIndex
int positionToIndex(double pos)
Convert image position to nearest integer index.
Definition: ImageUtils.h:69
iter
int iter
Definition: GaussianCentroid.cc:73
lsst::utils::isnan
int isnan(T t)
Definition: ieee.h:110
lsst::afw::geom::Box2I
An integer coordinate rectangle.
Definition: Box.h:53
image
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
d
int d
Definition: KDTree.cc:89
sum
double sum
Definition: NaiveFlux.cc:137
template<typename ImageT >
std::pair< double, double > lsst::meas::algorithms::detail::getFixedMomentsFlux ( ImageT const &  image,
double  bkgd,
double  xcen,
double  ycen,
detail::SdssShapeImpl const &  shape_ 
)

Return the flux of an object, using the aperture described by the SdssShape object.

The SdssShape algorithm calculates an elliptical Gaussian fit to an object, so the "aperture" is an elliptical Gaussian

Returns
A std::pair of the flux and its error
Parameters
imagethe data to process
bkgdbackground level
xcenx-centre of object
yceny-centre of object
shape_The SdssShape of the object

Definition at line 698 of file SdssShape.cc.

704 {
705  detail::SdssShapeImpl shape = shape_; // we need a mutable copy
706 
707  afwGeom::BoxI const& bbox = set_amom_bbox(image.getWidth(), image.getHeight(), xcen, ycen,
708  shape.getIxx(), shape.getIxy(), shape.getIyy());
709 
710  boost::tuple<std::pair<bool, double>, double, double, double> weights =
711  getWeights(shape.getIxx(), shape.getIxy(), shape.getIyy());
712  double const NaN = std::numeric_limits<double>::quiet_NaN();
713  if (!weights.get<0>().first) {
714  return std::make_pair(NaN, NaN);
715  }
716 
717  double const w11 = weights.get<1>();
718  double const w12 = weights.get<2>();
719  double const w22 = weights.get<3>();
720  bool const interp = shouldInterp(shape.getIxx(), shape.getIyy(), weights.get<0>().second);
721 
722  double I0 = 0; // amplitude of Gaussian
723  calcmom<true>(ImageAdaptor<ImageT>().getImage(image), xcen, ycen, bbox, bkgd, interp, w11, w12, w22,
724  &I0, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
725  /*
726  * We have enerything we need, but it isn't quite packaged right; we need an initialised SdssShapeImpl
727  */
728  shape.setI0(I0);
729 
730  {
731  int ix = static_cast<int>(xcen);
732  int iy = static_cast<int>(ycen);
733  float bkgd_var =
734  ImageAdaptor<ImageT>().getVariance(image, ix, iy); // XXX Overestimate as it includes object
735 
736  detail::SdssShapeImpl::Matrix4 fisher = calc_fisher(shape, bkgd_var); // Fisher matrix
737 
738  shape.setCovar(fisher.inverse());
739  }
740 
741  double const scale = shape.getFluxScale();
742  return std::make_pair(shape.getI0()*scale, shape.getI0Err()*scale);
743 }
lsst::afw::geom::Box2I
An integer coordinate rectangle.
Definition: Box.h:53
image
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
std::pair<double, double> lsst::meas::algorithms::detail::rotate ( double  x,
double  y,
double  angle 
)

Definition at line 432 of file SincFlux.cc.

432  {
433  double c = ::cos(angle);
434  double s = ::sin(angle);
435  return std::pair<double, double>(x*c + y*s, -x*s + y*c);
436  }
y
int y
Definition: GaussianCentroid.cc:36
x
int x
Definition: GaussianCentroid.cc:36

Variable Documentation

int const lsst::meas::algorithms::detail::SDSS_SHAPE_MAX_ITER = 100

Definition at line 12 of file SdssShape.h.

float const lsst::meas::algorithms::detail::SDSS_SHAPE_TOL1 = 1.0e-5

Definition at line 13 of file SdssShape.h.

float const lsst::meas::algorithms::detail::SDSS_SHAPE_TOL2 = 1.0e-4

Definition at line 14 of file SdssShape.h.

Generated on Mon Dec 8 2014 18:37:19 for LSSTApplications by   doxygen 1.8.5