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
SincCoeffs.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008-2013 LSST Corporation.
5  *
6  * This product includes software developed by the
7  * LSST Project (http://www.lsst.org/).
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the LSST License Statement and
20  * the GNU General Public License along with this program. If not,
21  * see <http://www.lsstcorp.org/LegalNotices/>.
22  */
23 
24 #include <complex>
25 
26 #include "boost/math/special_functions/bessel.hpp"
27 #include "boost/shared_array.hpp"
28 #include "fftw3.h"
29 
31 #include "lsst/afw/image/Image.h"
33 
34 namespace lsst { namespace meas { namespace base { namespace {
35 
36 // Convenient wrapper for a Bessel function
37 inline double J1(double const x) {
38  return boost::math::cyl_bessel_j(1, x);
39 }
40 
41 // sinc function
42 template<typename T>
43 inline T sinc(T const x) {
44  return (x != 0.0) ? (std::sin(x) / x) : 1.0;
45 }
46 
47 /*
48  * Define a circular aperture function object g_i, cos-tapered?
49  */
50 template<typename CoordT>
51 class CircularAperture {
52 public:
53 
54  CircularAperture(
55  CoordT const radius1,
56  CoordT const radius2,
57  CoordT const taperwidth
58  ):
59  _radius1(radius1),
60  _radius2(radius2),
61  _taperwidth1(taperwidth),
62  _taperwidth2(taperwidth),
63  _k1(1.0/(2.0*taperwidth)),
64  _k2(1.0/(2.0*taperwidth)),
65  _taperLo1(radius1 - 0.5*taperwidth),
66  _taperHi1(radius1 + 0.5*taperwidth),
67  _taperLo2(radius2 - 0.5*taperwidth),
68  _taperHi2(radius2 + 0.5*taperwidth) {
69 
70  // if we're asked for a radius smaller than our taperwidth,
71  // adjust the taper width smaller so it fits exactly
72  // with smooth derivative=0 at r=0
73 
74  if (_radius1 > _radius2) {
75  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
76  (boost::format("rad2 less than rad1: (rad1=%.2f, rad2=%.2f) ") %
77  _radius1 % _radius2).str());
78  }
79  if (_radius1 < 0.0 || _radius2 < 0.0) {
80  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
81  (boost::format("radii must be >= 0 (rad1=%.2f, rad2=%.2f) ") %
82  _radius1 % _radius2).str());
83  }
84 
85  if (_radius1 == 0) {
86  _taperwidth1 = 0.0;
87  _k1 = 0.0;
88  }
89 
90  // if we don't have room to taper at r=0
91  if ( _radius1 < 0.5*_taperwidth1) {
92  _taperwidth1 = 2.0*_radius1;
93  _k1 = 1.0/(2.0*_taperwidth1);
94  }
95 
96  // if we don't have room to taper between r1 and r2
97  if ((_radius2 - _radius1) < 0.5*(_taperwidth1+_taperwidth2)) {
98 
99  // if we *really* don't have room ... taper1 by itself is too big
100  // - set taper1,2 to be equal and split the r2-r1 range
101  if ((_radius2 - _radius2) < 0.5*_taperwidth1) {
103  _k1 = _k2 = 1.0/(2.0*_taperwidth1);
104 
105  // if there's room for taper1, but not taper1 and 2
106  } else {
108  _k2 = 1.0/(2.0*_taperwidth2);
109  }
110 
115  }
116  }
117 
118  // When called, return the throughput at the requested x,y
119  // todo: replace the sinusoid taper with a band-limited
120  CoordT operator() (CoordT const x, CoordT const y) const {
121  CoordT const xyrad = std::sqrt(x*x + y*y);
122  if ( xyrad < _taperLo1 ) {
123  return 0.0;
124  } else if (xyrad >= _taperLo1 && xyrad <= _taperHi1 ) {
125  return 0.5*(1.0 + std::cos( (afw::geom::TWOPI*_k1)*(xyrad - _taperHi1)));
126  } else if (xyrad > _taperHi1 && xyrad <= _taperLo2 ) {
127  return 1.0;
128  } else if (xyrad > _taperLo2 && xyrad <= _taperHi2 ) {
129  return 0.5*(1.0 + std::cos( (afw::geom::TWOPI*_k2)*(xyrad - _taperLo2)));
130  } else {
131  return 0.0;
132  }
133  }
134 
135  CoordT getRadius1() { return _radius1; }
136  CoordT getRadius2() { return _radius2; }
137 
138 private:
141  CoordT _k1, _k2; // the angular wavenumber corresponding to a cosine with wavelength 2*taperwidth
144 };
145 
146 
147 template<typename CoordT>
148 class CircApPolar : public std::unary_function<CoordT, CoordT> {
149 public:
150  CircApPolar(double radius, double taperwidth) : _ap(CircularAperture<CoordT>(0.0, radius, taperwidth)) {}
151  CoordT operator() (double r) const {
152  return r*_ap(r, 0.0);
153  }
154 private:
155  CircularAperture<CoordT> _ap;
156 };
157 
158 /*
159  * Define a Sinc functor to be integrated over for Sinc interpolation
160  */
161 template<typename IntegrandT>
162 class SincAperture : public std::binary_function<IntegrandT, IntegrandT, IntegrandT> {
163 public:
164 
165  SincAperture(
166  CircularAperture<IntegrandT> const &ap,
167  int const ix, // sinc center x
168  int const iy // sinc center y
169  )
170  : _ap(ap), _ix(ix), _iy(iy) {}
171 
172  IntegrandT operator() (IntegrandT const x, IntegrandT const y) const {
173  double const fourierConvention = afw::geom::PI;
174  double const dx = fourierConvention*(x - _ix);
175  double const dy = fourierConvention*(y - _iy);
176  double const fx = sinc<double>(dx);
177  double const fy = sinc<double>(dy);
178  return (1.0 + _ap(x, y)*fx*fy);
179  }
180 
181 private:
182  CircularAperture<IntegrandT> const &_ap;
183  double _ix, _iy;
184 };
185 
186 class GaussPowerFunctor : public std::binary_function<double, double, double> {
187 public:
188  GaussPowerFunctor(double sigma) : _sigma(sigma) {}
189 
190  double operator()(double kx, double ky) const {
191  double const k = ::sqrt(kx*kx + ky*ky);
192  double const gauss = ::exp(-0.5*k*k*_sigma*_sigma);
193  return gauss*gauss;
194  }
195 private:
196  double _sigma;
197 };
198 
199 std::pair<double, double> computeGaussLeakage(double const sigma) {
200 
201  GaussPowerFunctor gaussPower(sigma);
202 
203  double lim = afw::geom::PI;
204 
205  // total power: integrate GaussPowerFunctor -inf<x<inf, -inf<y<inf (can be done analytically)
206  double powerInf = afw::geom::PI/(sigma*sigma);
207 
208  // true power: integrate GaussPowerFunctor -lim<x<lim, -lim<y<lim (must be done numerically)
209  double truePower = afw::math::integrate2d(gaussPower, -lim, lim, -lim, lim, 1.0e-8);
210  double trueLeak = (powerInf - truePower)/powerInf;
211 
212  // estimated power: function is circular, but coords are cartesian
213  // - true power does the actual integral numerically, but we can estimate it by integrating
214  // in polar coords over lim <= radius < infinity. The integral is analytic.
215  double estLeak = ::exp(-sigma*sigma*afw::geom::PI*afw::geom::PI)/powerInf;
216 
217  return std::pair<double, double>(trueLeak, estLeak);
218 
219 }
220 
221 template<typename PixelT>
223 calcImageRealSpace(double const rad1, double const rad2, double const taperwidth) {
224 
225  PixelT initweight = 0.0; // initialize the coeff values
226 
227  int log2 = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
228  if (log2 < 3) { log2 = 3; }
229  int hwid = pow(2, log2);
230  int width = 2*hwid - 1;
231 
232  int const xwidth = width;
233  int const ywidth = width;
234 
235  int const x0 = -xwidth/2;
236  int const y0 = -ywidth/2;
237 
238  // create an image to hold the coefficient image
239  typename afw::image::Image<PixelT>::Ptr coeffImage =
240  boost::make_shared<afw::image::Image<PixelT> >(afw::geom::ExtentI(xwidth, ywidth), initweight);
241  coeffImage->setXY0(x0, y0);
242 
243  // create the aperture function object
244  // determine the radius to use that makes 'radius' the effective radius of the aperture
245  double tolerance = 1.0e-12;
246  double dr = 1.0e-6;
247  double err = 2.0*tolerance;
248  double apEff = afw::geom::PI*rad2*rad2;
249  double radIn = rad2;
250  int maxIt = 20;
251  int i = 0;
252  while (err > tolerance && i < maxIt) {
253  CircApPolar<double> apPolar1(radIn, taperwidth);
254  CircApPolar<double> apPolar2(radIn+dr, taperwidth);
255  double a1 = afw::geom::TWOPI * afw::math::integrate(apPolar1, 0.0, radIn+taperwidth, tolerance);
256  double a2 = afw::geom::TWOPI * afw::math::integrate(apPolar2, 0.0, radIn+dr+taperwidth, tolerance);
257  double dadr = (a2 - a1)/dr;
258  double radNew = radIn - (a1 - apEff)/dadr;
259  err = (a1 - apEff)/apEff;
260  radIn = radNew;
261  i++;
262  }
263  CircularAperture<double> ap(rad1, rad2, taperwidth);
264 
265  /* ******************************* */
266  // integrate over the aperture
267 
268  // the limits of the integration over the sinc aperture
269  double const limit = rad2 + taperwidth;
270  double const x1 = -limit;
271  double const x2 = limit;
272  double const y1 = -limit;
273  double const y2 = limit;
274 
275  for (int iY = y0; iY != y0 + coeffImage->getHeight(); ++iY) {
276  int iX = x0;
277  typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY-y0);
278  for (
279  typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY-y0);
280  ptr != end;
281  ++ptr
282  ) {
283 
284  // create a sinc function in the CircularAperture at our location
285  SincAperture<double> sincAp(ap, iX, iY);
286 
287  // integrate the sinc
288  PixelT integral = afw::math::integrate2d(sincAp, x1, x2, y1, y2, 1.0e-8);
289 
290  // we actually integrated function+1.0 and now must subtract the excess volume
291  // - just force it to zero in the corners
292  double const dx = iX;
293  double const dy = iY;
294  *ptr = (std::sqrt(dx*dx + dy*dy) < xwidth/2) ?
295  integral - (x2 - x1)*(y2 - y1) : 0.0;
296  ++iX;
297  }
298  }
299 
300 
301  double sum = 0.0;
302  for (int iY = y0; iY != y0 + coeffImage->getHeight(); ++iY) {
303  typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY-y0);
304  for (
305  typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY-y0);
306  ptr != end; ++ptr
307  ) {
308  sum += *ptr;
309  }
310  }
311 
312 #if 0 // debugging
313  coeffImage->writeFits("cimage.fits");
314 #endif
315 
316  return coeffImage;
317 }
318 
319 class FftShifter {
320 public:
321  FftShifter(int xwid) : _xwid(xwid) {}
322  int shift(int x) {
323  if (x >= _xwid/2) {
324  return x - _xwid/2;
325  } else {
326  return x + _xwid/2 + 1;
327  }
328  }
329 private:
330  int _xwid;
331 };
332 
333 std::pair<double, double> rotate(double x, double y, double angle) {
334  double c = ::cos(angle);
335  double s = ::sin(angle);
336  return std::pair<double, double>(x*c + y*s, -x*s + y*c);
337 }
338 
339 /* todo
340  * - try sub pixel shift if it doesn't break even symmetry
341  * - put values directly in an Image
342  * - precompute the plan
343  */
344 
345 template<typename PixelT>
346 typename afw::image::Image<PixelT>::Ptr calcImageKSpaceCplx(double const rad1, double const rad2,
347  double const posAng, double const ellipticity
348  ) {
349 
350  // we only need a half-width due to symmetry
351  // make the hwid 2*rad2 so we have some buffer space and round up to the next power of 2
352  int log2 = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
353  if (log2 < 3) { log2 = 3; }
354  int hwid = pow(2, log2);
355  int wid = 2*hwid - 1;
356  int xcen = wid/2, ycen = wid/2;
357  FftShifter fftshift(wid);
358 
359  boost::shared_array<std::complex<double> > cimg(new std::complex<double>[wid*wid]);
360  std::complex<double> *c = cimg.get();
361  // fftplan args: nx, ny, *in, *out, direction, flags
362  // - done in-situ if *in == *out
363  fftw_plan plan = fftw_plan_dft_2d(wid, wid,
364  reinterpret_cast<fftw_complex*>(c),
365  reinterpret_cast<fftw_complex*>(c),
366  FFTW_BACKWARD, FFTW_ESTIMATE);
367 
368  // compute the k-space values and put them in the cimg array
369  double const twoPiRad1 = afw::geom::TWOPI*rad1;
370  double const twoPiRad2 = afw::geom::TWOPI*rad2;
371  double const scale = (1.0 - ellipticity);
372  for (int iY = 0; iY < wid; ++iY) {
373  int const fY = fftshift.shift(iY);
374  double const ky = (static_cast<double>(iY) - ycen)/wid;
375 
376  for (int iX = 0; iX < wid; ++iX) {
377 
378  int const fX = fftshift.shift(iX);
379  double const kx = static_cast<double>(iX - xcen)/wid;
380 
381  // rotate
382  std::pair<double, double> coo = rotate(kx, ky, posAng);
383  double kxr = coo.first;
384  double kyr = coo.second;
385  // rescale
386  double const k = ::sqrt(kxr*kxr + scale*scale*kyr*kyr);
387 
388  double const airy1 = (rad1 > 0 ? rad1*J1(twoPiRad1*k) : 0.0)/k;
389  double const airy2 = rad2*J1(twoPiRad2*k)/k;
390  double const airy = airy2 - airy1;
391 
392  c[fY*wid + fX] = std::complex<double>(scale*airy, 0.0);
393  }
394  }
395  c[0] = scale*afw::geom::PI*(rad2*rad2 - rad1*rad1);
396 
397  // perform the fft and clean up after ourselves
398  fftw_execute(plan);
399  fftw_destroy_plan(plan);
400 
401  // put the coefficients into an image
402  typename afw::image::Image<PixelT>::Ptr coeffImage =
403  boost::make_shared<afw::image::Image<PixelT> >(afw::geom::ExtentI(wid, wid), 0.0);
404 
405  for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
406  int iX = 0;
407  typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
408  for (
409  typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY);\
410  ptr != end;
411  ++ptr
412  ) {
413  int fX = fftshift.shift(iX);
414  int fY = fftshift.shift(iY);
415  *ptr = static_cast<PixelT>(c[fY*wid + fX].real()/(wid*wid));
416  iX++;
417  }
418  }
419 
420  // reset the origin to be the middle of the image
421  coeffImage->setXY0(-wid/2, -wid/2);
422  return coeffImage;
423 }
424 
425 
426 
427 // I'm not sure why this doesn't work with DCT-I (REDFT00), the DCT should take advantage of symmetry
428 // and be much faster than the DFT. It runs but the numbers are slightly off ...
429 // but I have to do it as real-to-halfcplx (R2HC) to get the correct numbers.
430 
431 template<typename PixelT>
432 typename afw::image::Image<PixelT>::Ptr calcImageKSpaceReal(double const rad1, double const rad2) {
433 
434  // we only need a half-width due to symmertry
435  // make the hwid 2*rad2 so we have some buffer space and round up to the next power of 2
436  int log2 = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
437  if (log2 < 3) { log2 = 3; }
438  int hwid = pow(2, log2);
439  int wid = 2*hwid - 1;
440  int xcen = wid/2, ycen = wid/2;
441  FftShifter fftshift(wid);
442 
443  boost::shared_array<double> cimg(new double[wid*wid]);
444  double *c = cimg.get();
445  // fftplan args: nx, ny, *in, *out, kindx, kindy, flags
446  // - done in-situ if *in == *out
447  fftw_plan plan = fftw_plan_r2r_2d(wid, wid, c, c, FFTW_R2HC, FFTW_R2HC, FFTW_ESTIMATE);
448 
449  // compute the k-space values and put them in the cimg array
450  double const twoPiRad1 = afw::geom::TWOPI*rad1;
451  double const twoPiRad2 = afw::geom::TWOPI*rad2;
452  for (int iY = 0; iY < wid; ++iY) {
453 
454  int const fY = fftshift.shift(iY);
455  double const ky = (static_cast<double>(iY) - ycen)/wid;
456 
457  for (int iX = 0; iX < wid; ++iX) {
458  int const fX = fftshift.shift(iX);
459 
460  // emacs indent breaks if this isn't separte
461  double const iXcen = static_cast<double>(iX - xcen);
462  double const kx = iXcen/wid;
463 
464  double const k = ::sqrt(kx*kx + ky*ky);
465  double const airy1 = (rad1 > 0 ? rad1*J1(twoPiRad1*k) : 0.0)/k;
466  double const airy2 = rad2*J1(twoPiRad2*k)/k;
467  double const airy = airy2 - airy1;
468  c[fY*wid + fX] = airy;
469 
470  }
471  }
472  int fxy = fftshift.shift(wid/2);
473  c[fxy*wid + fxy] = afw::geom::PI*(rad2*rad2 - rad1*rad1);
474 
475  // perform the fft and clean up after ourselves
476  fftw_execute(plan);
477  fftw_destroy_plan(plan);
478 
479  // put the coefficients into an image
480  typename afw::image::Image<PixelT>::Ptr coeffImage =
481  boost::make_shared<afw::image::Image<PixelT> >(afw::geom::ExtentI(wid, wid), 0.0);
482 
483  for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
484  int iX = 0;
485  typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
486  for (
487  typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY);
488  ptr != end;
489  ++ptr
490  ) {
491 
492  // now need to reflect the quadrant we solved to the other three
493  int fX = iX < hwid ? hwid - iX - 1 : iX - hwid + 1;
494  int fY = iY < hwid ? hwid - iY - 1 : iY - hwid + 1;
495  *ptr = static_cast<PixelT>(c[fY*wid + fX]/(wid*wid));
496  iX++;
497  }
498  }
499 
500  // reset the origin to be the middle of the image
501  coeffImage->setXY0(-wid/2, -wid/2);
502  return coeffImage;
503 }
504 
505 } // anonymous
506 
507 template<typename PixelT>
509 {
510  static SincCoeffs<PixelT> instance;
511  return instance;
512 }
513 
514 template<typename PixelT>
515 void SincCoeffs<PixelT>::cache(float r1, float r2)
516 {
517  if (r1 < 0.0 || r2 < r1) {
518  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
519  (boost::format("Invalid r1,r2 = %f,%f") % r1 % r2).str());
520  }
521  double const innerFactor = r1/r2;
522  afw::geom::ellipses::Axes axes(r2, r2, 0.0);
523  if (!getInstance()._lookup(axes, innerFactor)) {
524  PTR(typename SincCoeffs<PixelT>::CoeffT) coeff = calculate(axes, innerFactor);
525  coeff->markPersistent();
526  getInstance()._cache[r2][innerFactor] = coeff;
527  }
528 }
529 
530 template<typename PixelT>
532 SincCoeffs<PixelT>::get(afw::geom::ellipses::Axes const& axes, float const innerFactor)
533 {
534  CONST_PTR(CoeffT) coeff = getInstance()._lookup(axes, innerFactor);
535  return coeff ? coeff : calculate(axes, innerFactor);
536 }
537 
538 template<typename PixelT>
540 SincCoeffs<PixelT>::_lookup(afw::geom::ellipses::Axes const& axes, double const innerFactor) const
541 {
542  if (innerFactor < 0.0 || innerFactor > 1.0) {
543  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
544  (boost::format("innerFactor = %f is not between 0 and 1") % innerFactor).str());
545  }
546 
548 
549  // We only cache circular apertures
550  if (!FuzzyCompare<float>().isEqual(axes.getA(), axes.getB())) {
551  return null;
552  }
553  typename CoeffMapMap::const_iterator iter1 = _cache.find(axes.getA());
554  if (iter1 == _cache.end()) {
555  return null;
556  }
557  typename CoeffMap::const_iterator iter2 = iter1->second.find(innerFactor);
558  return (iter2 == iter1->second.end()) ? null : iter2->second;
559 }
560 
561 template<typename PixelT>
563 SincCoeffs<PixelT>::calculate(afw::geom::ellipses::Axes const& axes, double const innerFactor)
564 {
565  if (innerFactor < 0.0 || innerFactor > 1.0) {
566  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
567  (boost::format("innerFactor = %f is not between 0 and 1") % innerFactor).str());
568  }
569 
570  // Kspace-real is fastest, but only slightly faster than kspace cplx
571  // but real won't work for elliptical apertures due to symmetries assumed for real transform
572 
573  double const rad1 = axes.getA() * innerFactor;
574  double const rad2 = axes.getA();
575  // if there's no angle and no ellipticity
576  if (FuzzyCompare<float>().isEqual(axes.getA(), axes.getB())) {
577  // here we call the real transform
578  return calcImageKSpaceReal<PixelT>(rad1, rad2);
579  } else {
580  // here we call the complex transform
581  double const ellipticity = 1.0 - axes.getB()/axes.getA();
582  return calcImageKSpaceCplx<PixelT>(rad1, rad2, axes.getTheta(), ellipticity);
583  }
584 }
585 
586 template class SincCoeffs<float>;
587 template class SincCoeffs<double>;
588 
589 }}} // namespace lsst::meas::base
590 
int y
CoordT _taperLo1
Definition: SincCoeffs.cc:142
#define PTR(...)
Definition: base.h:41
double _ix
Definition: SincCoeffs.cc:183
double _sigma
Definition: SincCoeffs.cc:196
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
#define CONST_PTR(...)
Definition: base.h:47
CoordT _radius2
Definition: SincCoeffs.cc:139
void shift(Vector< T, N > const &offset, Array< std::complex< T >, N, C > const &array, int const real_last_dim)
Perform a Fourier-space translation transform.
Definition: FourierOps.h:146
double const TWOPI
Definition: Angle.h:19
Extent< int, N > ceil(Extent< double, N > const &input)
boost::shared_ptr< Image< PixelT > > Ptr
Definition: Image.h:418
Extent< int, 2 > ExtentI
Definition: Extent.h:352
int const x0
Definition: saturated.cc:45
boost::enable_if< typename ExpressionTraits< Scalar >::IsScalar, Scalar >::type sum(Scalar const &scalar)
Definition: operators.h:1250
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
CoordT _taperLo2
Definition: SincCoeffs.cc:143
static SincCoeffs & getInstance()
Definition: SincCoeffs.cc:508
Support for 2-D images.
CircularAperture< CoordT > _ap
Definition: SincCoeffs.cc:155
double _iy
Definition: SincCoeffs.cc:183
static void cache(float rInner, float rOuter)
Definition: SincCoeffs.cc:515
int x
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a...
Definition: Axes.h:45
int _xwid
Definition: SincCoeffs.cc:330
CoordT _taperHi2
Definition: SincCoeffs.cc:143
_view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: Image.h:151
CoordT _radius1
Definition: SincCoeffs.cc:139
CoordT _taperwidth2
Definition: SincCoeffs.cc:140
double const PI
The ratio of a circle&#39;s circumference to diameter.
Definition: Angle.h:18
CoordT _taperHi1
Definition: SincCoeffs.cc:142
int const y0
Definition: saturated.cc:45
Compute 1d and 2d integral.
CoordT _taperwidth1
Definition: SincCoeffs.cc:140
CoordT _k2
Definition: SincCoeffs.cc:141
CoordT _k1
Definition: SincCoeffs.cc:141
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