LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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/geom/Angle.h"
32#include "lsst/geom/Extent.h"
35
36namespace lsst {
37namespace meas {
38namespace base {
39namespace {
40
41// Convenient wrapper for a Bessel function
42inline double J1(double const x) { return boost::math::cyl_bessel_j(1, x); }
43
44// sinc function
45template <typename T>
46inline T sinc(T const x) {
47 return (x != 0.0) ? (std::sin(x) / x) : 1.0;
48}
49
50/*
51 * Define a circular aperture function object g_i, cos-tapered?
52 */
53template <typename CoordT>
54class CircularAperture {
55public:
56 CircularAperture(
57 CoordT const radius1,
58 CoordT const radius2,
59 CoordT const taperwidth
60 )
61 : _radius1(radius1),
62 _radius2(radius2),
63 _taperwidth1(taperwidth),
64 _taperwidth2(taperwidth),
65 _k1(1.0 / (2.0 * taperwidth)),
66 _k2(1.0 / (2.0 * taperwidth)),
67 _taperLo1(radius1 - 0.5 * taperwidth),
68 _taperHi1(radius1 + 0.5 * taperwidth),
69 _taperLo2(radius2 - 0.5 * taperwidth),
70 _taperHi2(radius2 + 0.5 * taperwidth) {
71 // if we're asked for a radius smaller than our taperwidth,
72 // adjust the taper width smaller so it fits exactly
73 // with smooth derivative=0 at r=0
74
75 if (_radius1 > _radius2) {
76 throw LSST_EXCEPT(
77 pex::exceptions::InvalidParameterError,
78 (boost::format("rad2 less than rad1: (rad1=%.2f, rad2=%.2f) ") % _radius1 % _radius2)
79 .str());
80 }
81 if (_radius1 < 0.0 || _radius2 < 0.0) {
82 throw LSST_EXCEPT(
83 pex::exceptions::InvalidParameterError,
84 (boost::format("radii must be >= 0 (rad1=%.2f, rad2=%.2f) ") % _radius1 % _radius2)
85 .str());
86 }
87
88 if (_radius1 == 0) {
89 _taperwidth1 = 0.0;
90 _k1 = 0.0;
91 }
92
93 // if we don't have room to taper at r=0
94 if (_radius1 < 0.5 * _taperwidth1) {
95 _taperwidth1 = 2.0 * _radius1;
96 _k1 = 1.0 / (2.0 * _taperwidth1);
97 }
98
99 // if we don't have room to taper between r1 and r2
100 if ((_radius2 - _radius1) < 0.5 * (_taperwidth1 + _taperwidth2)) {
101 // if we *really* don't have room ... taper1 by itself is too big
102 // - set taper1,2 to be equal and split the r2-r1 range
103 if ((_radius2 - _radius2) < 0.5 * _taperwidth1) {
104 _taperwidth1 = _taperwidth2 = 0.5 * (_radius2 - _radius1);
105 _k1 = _k2 = 1.0 / (2.0 * _taperwidth1);
106
107 // if there's room for taper1, but not taper1 and 2
108 } else {
109 _taperwidth2 = _radius2 - _radius1 - _taperwidth1;
110 _k2 = 1.0 / (2.0 * _taperwidth2);
111 }
112
113 _taperLo1 = _radius1 - 0.5 * _taperwidth1;
114 _taperHi1 = _radius1 + 0.5 * _taperwidth1;
115 _taperLo2 = _radius2 - 0.5 * _taperwidth2;
116 _taperHi2 = _radius2 + 0.5 * _taperwidth2;
117 }
118 }
119
120 // When called, return the throughput at the requested x,y
121 // todo: replace the sinusoid taper with a band-limited
122 CoordT operator()(CoordT const x, CoordT const y) const {
123 CoordT const xyrad = std::sqrt(x * x + y * y);
124 if (xyrad < _taperLo1) {
125 return 0.0;
126 } else if (xyrad >= _taperLo1 && xyrad <= _taperHi1) {
127 return 0.5 * (1.0 + std::cos((geom::TWOPI * _k1) * (xyrad - _taperHi1)));
128 } else if (xyrad > _taperHi1 && xyrad <= _taperLo2) {
129 return 1.0;
130 } else if (xyrad > _taperLo2 && xyrad <= _taperHi2) {
131 return 0.5 * (1.0 + std::cos((geom::TWOPI * _k2) * (xyrad - _taperLo2)));
132 } else {
133 return 0.0;
134 }
135 }
136
137 CoordT getRadius1() { return _radius1; }
138 CoordT getRadius2() { return _radius2; }
139
140private:
141 CoordT _radius1, _radius2;
142 CoordT _taperwidth1, _taperwidth2;
143 CoordT _k1, _k2; // the angular wavenumber corresponding to a cosine with wavelength 2*taperwidth
144 CoordT _taperLo1, _taperHi1;
145 CoordT _taperLo2, _taperHi2;
146};
147
148template <typename CoordT>
149class CircApPolar {
150public:
151 CircApPolar(double radius, double taperwidth) : _ap(CircularAperture<CoordT>(0.0, radius, taperwidth)) {}
152 CoordT operator()(double r) const { return r * _ap(r, 0.0); }
153
154private:
155 CircularAperture<CoordT> _ap;
156};
157
158/*
159 * Define a Sinc functor to be integrated over for Sinc interpolation
160 */
161template <typename IntegrandT>
162class SincAperture {
163public:
164 SincAperture(CircularAperture<IntegrandT> const& ap,
165 int const ix, // sinc center x
166 int const iy // sinc center y
167 )
168 : _ap(ap), _ix(ix), _iy(iy) {}
169
170 IntegrandT operator()(IntegrandT const x, IntegrandT const y) const {
171 double const fourierConvention = geom::PI;
172 double const dx = fourierConvention * (x - _ix);
173 double const dy = fourierConvention * (y - _iy);
174 double const fx = sinc<double>(dx);
175 double const fy = sinc<double>(dy);
176 return (1.0 + _ap(x, y) * fx * fy);
177 }
178
179private:
180 CircularAperture<IntegrandT> const& _ap;
181 double _ix, _iy;
182};
183
184class GaussPowerFunctor {
185public:
186 GaussPowerFunctor(double sigma) : _sigma(sigma) {}
187
188 double operator()(double kx, double ky) const {
189 double const k = ::sqrt(kx * kx + ky * ky);
190 double const gauss = ::exp(-0.5 * k * k * _sigma * _sigma);
191 return gauss * gauss;
192 }
193
194private:
195 double _sigma;
196};
197
198std::pair<double, double> computeGaussLeakage(double const sigma) {
199 GaussPowerFunctor gaussPower(sigma);
200
201 double lim = geom::PI;
202
203 // total power: integrate GaussPowerFunctor -inf<x<inf, -inf<y<inf (can be done analytically)
204 double powerInf = geom::PI / (sigma * sigma);
205
206 // true power: integrate GaussPowerFunctor -lim<x<lim, -lim<y<lim (must be done numerically)
207 double truePower = afw::math::integrate2d(gaussPower, -lim, lim, -lim, lim, 1.0e-8);
208 double trueLeak = (powerInf - truePower) / powerInf;
209
210 // estimated power: function is circular, but coords are cartesian
211 // - true power does the actual integral numerically, but we can estimate it by integrating
212 // in polar coords over lim <= radius < infinity. The integral is analytic.
213 double estLeak = ::exp(-sigma * sigma * geom::PI * geom::PI) / powerInf;
214
215 return std::pair<double, double>(trueLeak, estLeak);
216}
217
218template <typename PixelT>
219std::shared_ptr<afw::image::Image<PixelT>> calcImageRealSpace(double const rad1, double const rad2,
220 double const taperwidth) {
221 PixelT initweight = 0.0; // initialize the coeff values
222
223 int log2 = static_cast<int>(::ceil(::log10(2.0 * rad2) / log10(2.0)));
224 if (log2 < 3) {
225 log2 = 3;
226 }
227 int hwid = pow(2, log2);
228 int width = 2 * hwid - 1;
229
230 int const xwidth = width;
231 int const ywidth = width;
232
233 int const x0 = -xwidth / 2;
234 int const y0 = -ywidth / 2;
235
236 // create an image to hold the coefficient image
237 auto coeffImage = std::make_shared<afw::image::Image<PixelT>>(geom::ExtentI(xwidth, ywidth), initweight);
238 coeffImage->setXY0(x0, y0);
239
240 // create the aperture function object
241 // determine the radius to use that makes 'radius' the effective radius of the aperture
242 double tolerance = 1.0e-12;
243 double dr = 1.0e-6;
244 double err = 2.0 * tolerance;
245 double apEff = geom::PI * rad2 * rad2;
246 double radIn = rad2;
247 int maxIt = 20;
248 int i = 0;
249 while (err > tolerance && i < maxIt) {
250 CircApPolar<double> apPolar1(radIn, taperwidth);
251 CircApPolar<double> apPolar2(radIn + dr, taperwidth);
252 double a1 = geom::TWOPI * afw::math::integrate(apPolar1, 0.0, radIn + taperwidth, tolerance);
253 double a2 = geom::TWOPI * afw::math::integrate(apPolar2, 0.0, radIn + dr + taperwidth, tolerance);
254 double dadr = (a2 - a1) / dr;
255 double radNew = radIn - (a1 - apEff) / dadr;
256 err = (a1 - apEff) / apEff;
257 radIn = radNew;
258 i++;
259 }
260 CircularAperture<double> ap(rad1, rad2, taperwidth);
261
262 /* ******************************* */
263 // integrate over the aperture
264
265 // the limits of the integration over the sinc aperture
266 double const limit = rad2 + taperwidth;
267 double const x1 = -limit;
268 double const x2 = limit;
269 double const y1 = -limit;
270 double const y2 = limit;
271
272 for (int iY = y0; iY != y0 + coeffImage->getHeight(); ++iY) {
273 int iX = x0;
274 typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY - y0);
275 for (typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY - y0); ptr != end;
276 ++ptr) {
277 // create a sinc function in the CircularAperture at our location
278 SincAperture<double> sincAp(ap, iX, iY);
279
280 // integrate the sinc
281 PixelT integral = afw::math::integrate2d(sincAp, x1, x2, y1, y2, 1.0e-8);
282
283 // we actually integrated function+1.0 and now must subtract the excess volume
284 // - just force it to zero in the corners
285 double const dx = iX;
286 double const dy = iY;
287 *ptr = (std::sqrt(dx * dx + dy * dy) < xwidth / 2) ? integral - (x2 - x1) * (y2 - y1) : 0.0;
288 ++iX;
289 }
290 }
291
292 double sum = 0.0;
293 for (int iY = y0; iY != y0 + coeffImage->getHeight(); ++iY) {
294 typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY - y0);
295 for (typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY - y0); ptr != end;
296 ++ptr) {
297 sum += *ptr;
298 }
299 }
300
301#if 0 // debugging
302 coeffImage->writeFits("cimage.fits");
303#endif
304
305 return coeffImage;
306}
307
308class FftShifter {
309public:
310 FftShifter(int xwid) : _xwid(xwid) {}
311 int shift(int x) {
312 if (x >= _xwid / 2) {
313 return x - _xwid / 2;
314 } else {
315 return x + _xwid / 2 + 1;
316 }
317 }
318
319private:
320 int _xwid;
321};
322
323std::pair<double, double> rotate(double x, double y, double angle) {
324 double c = ::cos(angle);
325 double s = ::sin(angle);
326 return std::pair<double, double>(x * c + y * s, -x * s + y * c);
327}
328
329/* todo
330 * - try sub pixel shift if it doesn't break even symmetry
331 * - put values directly in an Image
332 * - precompute the plan
333 */
334
335template <typename PixelT>
336std::shared_ptr<afw::image::Image<PixelT>> calcImageKSpaceCplx(double const rad1, double const rad2,
337 double const posAng,
338 double const ellipticity) {
339 // we only need a half-width due to symmetry
340 // make the hwid 2*rad2 so we have some buffer space and round up to the next power of 2
341 int log2 = static_cast<int>(::ceil(::log10(2.0 * rad2) / log10(2.0)));
342 if (log2 < 3) {
343 log2 = 3;
344 }
345 int hwid = pow(2, log2);
346 int wid = 2 * hwid - 1;
347 int xcen = wid / 2, ycen = wid / 2;
348 FftShifter fftshift(wid);
349
350 boost::shared_array<std::complex<double>> cimg(new std::complex<double>[wid * wid]);
351 std::complex<double>* c = cimg.get();
352 // fftplan args: nx, ny, *in, *out, direction, flags
353 // - done in-situ if *in == *out
354 fftw_plan plan = fftw_plan_dft_2d(wid, wid, reinterpret_cast<fftw_complex*>(c),
355 reinterpret_cast<fftw_complex*>(c), FFTW_BACKWARD, FFTW_ESTIMATE);
356
357 // compute the k-space values and put them in the cimg array
358 double const twoPiRad1 = geom::TWOPI * rad1;
359 double const twoPiRad2 = geom::TWOPI * rad2;
360 double const scale = (1.0 - ellipticity);
361 for (int iY = 0; iY < wid; ++iY) {
362 int const fY = fftshift.shift(iY);
363 double const ky = (static_cast<double>(iY) - ycen) / wid;
364
365 for (int iX = 0; iX < wid; ++iX) {
366 int const fX = fftshift.shift(iX);
367 double const kx = static_cast<double>(iX - xcen) / wid;
368
369 // rotate
370 std::pair<double, double> coo = rotate(kx, ky, posAng);
371 double kxr = coo.first;
372 double kyr = coo.second;
373 // rescale
374 double const k = ::sqrt(kxr * kxr + scale * scale * kyr * kyr);
375
376 double const airy1 = (rad1 > 0 ? rad1 * J1(twoPiRad1 * k) : 0.0) / k;
377 double const airy2 = rad2 * J1(twoPiRad2 * k) / k;
378 double const airy = airy2 - airy1;
379
380 c[fY * wid + fX] = std::complex<double>(scale * airy, 0.0);
381 }
382 }
383 c[0] = scale * geom::PI * (rad2 * rad2 - rad1 * rad1);
384
385 // perform the fft and clean up after ourselves
386 fftw_execute(plan);
387 fftw_destroy_plan(plan);
388
389 // put the coefficients into an image
390 auto coeffImage = std::make_shared<afw::image::Image<PixelT>>(geom::ExtentI(wid, wid), 0.0);
391
392 for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
393 int iX = 0;
394 typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
395 for (typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY); ptr != end;
396 ++ptr) {
397 int fX = fftshift.shift(iX);
398 int fY = fftshift.shift(iY);
399 *ptr = static_cast<PixelT>(c[fY * wid + fX].real() / (wid * wid));
400 iX++;
401 }
402 }
403
404 // reset the origin to be the middle of the image
405 coeffImage->setXY0(-wid / 2, -wid / 2);
406 return coeffImage;
407}
408
409// I'm not sure why this doesn't work with DCT-I (REDFT00), the DCT should take advantage of symmetry
410// and be much faster than the DFT. It runs but the numbers are slightly off ...
411// but I have to do it as real-to-halfcplx (R2HC) to get the correct numbers.
412
413template <typename PixelT>
414std::shared_ptr<afw::image::Image<PixelT>> calcImageKSpaceReal(double const rad1, double const rad2) {
415 // we only need a half-width due to symmertry
416 // make the hwid 2*rad2 so we have some buffer space and round up to the next power of 2
417 int log2 = static_cast<int>(::ceil(::log10(2.0 * rad2) / log10(2.0)));
418 if (log2 < 3) {
419 log2 = 3;
420 }
421 int hwid = pow(2, log2);
422 int wid = 2 * hwid - 1;
423 int xcen = wid / 2, ycen = wid / 2;
424 FftShifter fftshift(wid);
425
426 boost::shared_array<std::complex<double>> cimg(new std::complex<double>[wid * wid]);
427 std::complex<double>* c = cimg.get();
428 // fftplan args: nx, ny, *in, *out, kindx, kindy, flags
429 // - done in-situ if *in == *out
430
431 fftw_plan plan = fftw_plan_dft_2d(wid, wid, reinterpret_cast<fftw_complex*>(c),
432 reinterpret_cast<fftw_complex*>(c), FFTW_BACKWARD, FFTW_ESTIMATE);
433 // compute the k-space values and put them in the cimg array
434 double const twoPiRad1 = geom::TWOPI * rad1;
435 double const twoPiRad2 = geom::TWOPI * rad2;
436 for (int iY = 0; iY < wid; ++iY) {
437 int const fY = fftshift.shift(iY);
438 double const ky = (static_cast<double>(iY) - ycen) / wid;
439
440 for (int iX = 0; iX < wid; ++iX) {
441 int const fX = fftshift.shift(iX);
442
443 // emacs indent breaks if this isn't separte
444 double const iXcen = static_cast<double>(iX - xcen);
445 double const kx = iXcen / wid;
446
447 double const k = ::sqrt(kx * kx + ky * ky);
448 double const airy1 = (rad1 > 0 ? rad1 * J1(twoPiRad1 * k) : 0.0) / k;
449 double const airy2 = rad2 * J1(twoPiRad2 * k) / k;
450 double const airy = airy2 - airy1;
451 c[fY * wid + fX] = airy;
452 }
453 }
454 int fxy = fftshift.shift(wid / 2);
455 c[fxy * wid + fxy] = {geom::PI * (rad2 * rad2 - rad1 * rad1), 0.};
456
457 // perform the fft and clean up after ourselves
458 fftw_execute(plan);
459 fftw_destroy_plan(plan);
460
461 // put the coefficients into an image
462 auto coeffImage = std::make_shared<afw::image::Image<PixelT>>(geom::ExtentI(wid, wid), 0.0);
463
464 for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
465 int iX = 0;
466 typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
467 for (typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY); ptr != end;
468 ++ptr) {
469 // now need to reflect the quadrant we solved to the other three
470 int fX = iX < hwid ? hwid - iX - 1 : iX - hwid + 1;
471 int fY = iY < hwid ? hwid - iY - 1 : iY - hwid + 1;
472 *ptr = static_cast<PixelT>(c[fY * wid + fX].real() / (wid * wid));
473 iX++;
474 }
475 }
476
477 // reset the origin to be the middle of the image
478 coeffImage->setXY0(-wid / 2, -wid / 2);
479 return coeffImage;
480}
481
482} // namespace
483
484template <typename PixelT>
485SincCoeffs<PixelT>& SincCoeffs<PixelT>::getInstance() {
486 static SincCoeffs<PixelT> instance;
487 return instance;
488}
489
490template <typename PixelT>
491void SincCoeffs<PixelT>::cache(float r1, float r2) {
492 if (r1 < 0.0 || r2 < r1) {
494 (boost::format("Invalid r1,r2 = %f,%f") % r1 % r2).str());
495 }
496 double const innerFactor = r1 / r2;
497 afw::geom::ellipses::Axes axes(r2, r2, 0.0);
498 if (!getInstance()._lookup(axes, innerFactor)) {
499 std::shared_ptr<typename SincCoeffs<PixelT>::CoeffT> coeff = calculate(axes, innerFactor);
500 getInstance()._cache[r2][innerFactor] = coeff;
501 }
502}
503
504template <typename PixelT>
506SincCoeffs<PixelT>::get(afw::geom::ellipses::Axes const& axes, float const innerFactor) {
507 std::shared_ptr<CoeffT const> coeff = getInstance()._lookup(axes, innerFactor);
508 return coeff ? coeff : calculate(axes, innerFactor);
509}
510
511template <typename PixelT>
513SincCoeffs<PixelT>::_lookup(afw::geom::ellipses::Axes const& axes, double const innerFactor) const {
514 if (innerFactor < 0.0 || innerFactor > 1.0) {
516 (boost::format("innerFactor = %f is not between 0 and 1") % innerFactor).str());
517 }
518
520
521 // We only cache circular apertures
522 if (!FuzzyCompare<float>().isEqual(axes.getA(), axes.getB())) {
523 return null;
524 }
525 typename CoeffMapMap::const_iterator iter1 = _cache.find(axes.getA());
526 if (iter1 == _cache.end()) {
527 return null;
528 }
529 typename CoeffMap::const_iterator iter2 = iter1->second.find(innerFactor);
530 return (iter2 == iter1->second.end()) ? null : iter2->second;
531}
532
533template <typename PixelT>
535SincCoeffs<PixelT>::calculate(afw::geom::ellipses::Axes const& axes, double const innerFactor) {
536 if (innerFactor < 0.0 || innerFactor > 1.0) {
538 (boost::format("innerFactor = %f is not between 0 and 1") % innerFactor).str());
539 }
540
541 // Kspace-real is fastest, but only slightly faster than kspace cplx
542 // but real won't work for elliptical apertures due to symmetries assumed for real transform
543
544 double const rad1 = axes.getA() * innerFactor;
545 double const rad2 = axes.getA();
546 // if there's no angle and no ellipticity
547 if (FuzzyCompare<float>().isEqual(axes.getA(), axes.getB())) {
548 // here we call the real transform
549 return calcImageKSpaceReal<PixelT>(rad1, rad2);
550 } else {
551 // here we call the complex transform
552 double const ellipticity = 1.0 - axes.getB() / axes.getA();
553 return calcImageKSpaceCplx<PixelT>(rad1, rad2, axes.getTheta(), ellipticity);
554 }
555}
556
557template class SincCoeffs<float>;
558template class SincCoeffs<double>;
559
560} // namespace base
561} // namespace meas
562} // namespace lsst
int end
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
table::Key< double > angle
afw::table::Key< double > sigma
uint64_t * ptr
Definition RangeSet.cc:95
int y
Definition SpanSet.cc:48
An ellipse core for the semimajor/semiminor axis and position angle parametrization (a,...
Definition Axes.h:47
double const getTheta() const
Definition Axes.h:57
double const getA() const
Definition Axes.h:51
double const getB() const
Definition Axes.h:54
typename _view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition ImageBase.h:133
A singleton to calculate and cache the coefficients for sinc photometry.
Definition SincCoeffs.h:45
static std::shared_ptr< CoeffT > calculate(afw::geom::ellipses::Axes const &outerEllipse, double const innerFactor=0.0)
Calculate the coefficients for an aperture.
static void cache(float rInner, float rOuter)
Cache the coefficients for a particular aperture.
static std::shared_ptr< CoeffT const > get(afw::geom::ellipses::Axes const &outerEllipse, float const innerRadiusFactor=0.0)
Get the coefficients for an aperture.
Reports invalid arguments.
Definition Runtime.h:66
T cos(T... args)
scale(algorithm, min, max=None, frame=None)
Definition ds9.py:108
auto integrate(UnaryFunctionT func, Arg const a, Arg const b, double eps=1.0e-6)
The 1D integrator.
Definition Integrate.h:765
auto integrate2d(BinaryFunctionT func, X x1, X x2, Y y1, Y y2, double eps=1.0e-6)
The 2D integrator.
Definition Integrate.h:779
double constexpr TWOPI
Definition Angle.h:41
Extent< int, 2 > ExtentI
Definition Extent.h:396
double constexpr PI
The ratio of a circle's circumference to diameter.
Definition Angle.h:40
uint8_t log2(uint64_t x)
Definition curve.h:105
T pow(T... args)
T real(T... args)
T rotate(T... args)
T sin(T... args)
T sqrt(T... args)
table::Key< table::Array< double > > coeff
Definition PsfexPsf.cc:362