LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
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"
33 #include "lsst/afw/image/Image.h"
35 
36 namespace lsst {
37 namespace meas {
38 namespace base {
39 namespace {
40 
41 // Convenient wrapper for a Bessel function
42 inline double J1(double const x) { return boost::math::cyl_bessel_j(1, x); }
43 
44 // sinc function
45 template <typename T>
46 inline 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  */
53 template <typename CoordT>
54 class CircularAperture {
55 public:
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 
140 private:
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 
148 template <typename CoordT>
149 class CircApPolar {
150 public:
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 
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 {
163 public:
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 
179 private:
180  CircularAperture<IntegrandT> const& _ap;
181  double _ix, _iy;
182 };
183 
184 class GaussPowerFunctor {
185 public:
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 
194 private:
195  double _sigma;
196 };
197 
198 std::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 
218 template <typename PixelT>
219 std::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 
308 class FftShifter {
309 public:
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 
319 private:
320  int _xwid;
321 };
322 
323 std::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 
335 template <typename PixelT>
336 std::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 
413 template <typename PixelT>
414 std::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<double> cimg(new double[wid * wid]);
427  double* c = cimg.get();
428  // fftplan args: nx, ny, *in, *out, kindx, kindy, flags
429  // - done in-situ if *in == *out
430  fftw_plan plan = fftw_plan_r2r_2d(wid, wid, c, c, FFTW_R2HC, FFTW_R2HC, FFTW_ESTIMATE);
431 
432  // compute the k-space values and put them in the cimg array
433  double const twoPiRad1 = geom::TWOPI * rad1;
434  double const twoPiRad2 = geom::TWOPI * rad2;
435  for (int iY = 0; iY < wid; ++iY) {
436  int const fY = fftshift.shift(iY);
437  double const ky = (static_cast<double>(iY) - ycen) / wid;
438 
439  for (int iX = 0; iX < wid; ++iX) {
440  int const fX = fftshift.shift(iX);
441 
442  // emacs indent breaks if this isn't separte
443  double const iXcen = static_cast<double>(iX - xcen);
444  double const kx = iXcen / wid;
445 
446  double const k = ::sqrt(kx * kx + ky * ky);
447  double const airy1 = (rad1 > 0 ? rad1 * J1(twoPiRad1 * k) : 0.0) / k;
448  double const airy2 = rad2 * J1(twoPiRad2 * k) / k;
449  double const airy = airy2 - airy1;
450  c[fY * wid + fX] = airy;
451  }
452  }
453  int fxy = fftshift.shift(wid / 2);
454  c[fxy * wid + fxy] = geom::PI * (rad2 * rad2 - rad1 * rad1);
455 
456  // perform the fft and clean up after ourselves
457  fftw_execute(plan);
458  fftw_destroy_plan(plan);
459 
460  // put the coefficients into an image
461  auto coeffImage = std::make_shared<afw::image::Image<PixelT>>(geom::ExtentI(wid, wid), 0.0);
462 
463  for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
464  int iX = 0;
465  typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
466  for (typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY); ptr != end;
467  ++ptr) {
468  // now need to reflect the quadrant we solved to the other three
469  int fX = iX < hwid ? hwid - iX - 1 : iX - hwid + 1;
470  int fY = iY < hwid ? hwid - iY - 1 : iY - hwid + 1;
471  *ptr = static_cast<PixelT>(c[fY * wid + fX] / (wid * wid));
472  iX++;
473  }
474  }
475 
476  // reset the origin to be the middle of the image
477  coeffImage->setXY0(-wid / 2, -wid / 2);
478  return coeffImage;
479 }
480 
481 } // namespace
482 
483 template <typename PixelT>
484 SincCoeffs<PixelT>& SincCoeffs<PixelT>::getInstance() {
485  static SincCoeffs<PixelT> instance;
486  return instance;
487 }
488 
489 template <typename PixelT>
490 void SincCoeffs<PixelT>::cache(float r1, float r2) {
491  if (r1 < 0.0 || r2 < r1) {
493  (boost::format("Invalid r1,r2 = %f,%f") % r1 % r2).str());
494  }
495  double const innerFactor = r1 / r2;
496  afw::geom::ellipses::Axes axes(r2, r2, 0.0);
497  if (!getInstance()._lookup(axes, innerFactor)) {
498  std::shared_ptr<typename SincCoeffs<PixelT>::CoeffT> coeff = calculate(axes, innerFactor);
499  getInstance()._cache[r2][innerFactor] = coeff;
500  }
501 }
502 
503 template <typename PixelT>
505 SincCoeffs<PixelT>::get(afw::geom::ellipses::Axes const& axes, float const innerFactor) {
506  std::shared_ptr<CoeffT const> coeff = getInstance()._lookup(axes, innerFactor);
507  return coeff ? coeff : calculate(axes, innerFactor);
508 }
509 
510 template <typename PixelT>
512 SincCoeffs<PixelT>::_lookup(afw::geom::ellipses::Axes const& axes, double const innerFactor) const {
513  if (innerFactor < 0.0 || innerFactor > 1.0) {
515  (boost::format("innerFactor = %f is not between 0 and 1") % innerFactor).str());
516  }
517 
519 
520  // We only cache circular apertures
521  if (!FuzzyCompare<float>().isEqual(axes.getA(), axes.getB())) {
522  return null;
523  }
524  typename CoeffMapMap::const_iterator iter1 = _cache.find(axes.getA());
525  if (iter1 == _cache.end()) {
526  return null;
527  }
528  typename CoeffMap::const_iterator iter2 = iter1->second.find(innerFactor);
529  return (iter2 == iter1->second.end()) ? null : iter2->second;
530 }
531 
532 template <typename PixelT>
534 SincCoeffs<PixelT>::calculate(afw::geom::ellipses::Axes const& axes, double const innerFactor) {
535  if (innerFactor < 0.0 || innerFactor > 1.0) {
537  (boost::format("innerFactor = %f is not between 0 and 1") % innerFactor).str());
538  }
539 
540  // Kspace-real is fastest, but only slightly faster than kspace cplx
541  // but real won't work for elliptical apertures due to symmetries assumed for real transform
542 
543  double const rad1 = axes.getA() * innerFactor;
544  double const rad2 = axes.getA();
545  // if there's no angle and no ellipticity
546  if (FuzzyCompare<float>().isEqual(axes.getA(), axes.getB())) {
547  // here we call the real transform
548  return calcImageKSpaceReal<PixelT>(rad1, rad2);
549  } else {
550  // here we call the complex transform
551  double const ellipticity = 1.0 - axes.getB() / axes.getA();
552  return calcImageKSpaceCplx<PixelT>(rad1, rad2, axes.getTheta(), ellipticity);
553  }
554 }
555 
556 template class SincCoeffs<float>;
557 template class SincCoeffs<double>;
558 
559 } // namespace base
560 } // namespace meas
561 } // namespace lsst
int end
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
table::Key< double > angle
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:49
uint64_t * ptr
Definition: RangeSet.cc:88
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.
Definition: SincCoeffs.cc:534
static void cache(float rInner, float rOuter)
Cache the coefficients for a particular aperture.
Definition: SincCoeffs.cc:490
static std::shared_ptr< CoeffT const > get(afw::geom::ellipses::Axes const &outerEllipse, float const innerRadiusFactor=0.0)
Get the coefficients for an aperture.
Definition: SincCoeffs.cc:505
Reports invalid arguments.
Definition: Runtime.h:66
T cos(T... args)
T log10(T... args)
def 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:767
auto integrate2d(BinaryFunctionT func, X x1, X x2, Y y1, Y y2, double eps=1.0e-6)
The 2D integrator.
Definition: Integrate.h:781
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
constexpr double TWOPI
Definition: Angle.h:40
Extent< int, 2 > ExtentI
Definition: Extent.h:396
Extent< int, N > ceil(Extent< double, N > const &input) noexcept
Return the component-wise ceil (round towards more positive).
Definition: Extent.cc:118
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
double sin(Angle const &a)
Definition: Angle.h:102
double cos(Angle const &a)
Definition: Angle.h:103
uint8_t log2(uint64_t x)
Definition: curve.h:98
A base class for image defects.
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