SincCoeffs.cc

Go to the documentation of this file.

// -*- lsst-c++ -*-
/*
 * LSST Data Management System
 * Copyright 2008-2013 LSST Corporation.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */

#include <complex>

#include "boost/math/special_functions/bessel.hpp"
#include "boost/shared_array.hpp"
#include "fftw3.h"

#include "lsst/meas/base/SincCoeffs.h"
#include "lsst/afw/image/Image.h"
#include "lsst/afw/math/Integrate.h"

namespace lsst { namespace meas { namespace base { namespace {

// Convenient wrapper for a Bessel function
inline double J1(double const x) {
    return boost::math::cyl_bessel_j(1, x);
}

// sinc function
template<typename T>
inline T sinc(T const x) {
    return (x != 0.0) ? (std::sin(x) / x) : 1.0;
}

/*
 * Define a circular aperture function object g_i, cos-tapered?
 */
template<typename CoordT>
class CircularAperture {
public:

    CircularAperture(
        CoordT const radius1,    
        CoordT const radius2,    
        CoordT const taperwidth 
    ):
        _radius1(radius1),
        _radius2(radius2),
        _taperwidth1(taperwidth),
        _taperwidth2(taperwidth),
        _k1(1.0/(2.0*taperwidth)),
        _k2(1.0/(2.0*taperwidth)),
        _taperLo1(radius1 - 0.5*taperwidth),
        _taperHi1(radius1 + 0.5*taperwidth),
        _taperLo2(radius2 - 0.5*taperwidth),
        _taperHi2(radius2 + 0.5*taperwidth) {

        // if we're asked for a radius smaller than our taperwidth,
        // adjust the taper width smaller so it fits exactly
        // with smooth derivative=0 at r=0

        if (_radius1 > _radius2) {
            throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                              (boost::format("rad2 less than rad1: (rad1=%.2f, rad2=%.2f) ") %
                               _radius1 % _radius2).str());
        }
        if (_radius1 < 0.0 || _radius2 < 0.0) {
            throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                              (boost::format("radii must be >= 0 (rad1=%.2f, rad2=%.2f) ") %
                               _radius1 % _radius2).str());
        }

        if (_radius1 == 0) {
            _taperwidth1 = 0.0;
            _k1 = 0.0;
        }

        // if we don't have room to taper at r=0
        if ( _radius1 < 0.5*_taperwidth1) {
            _taperwidth1 = 2.0*_radius1;
            _k1 = 1.0/(2.0*_taperwidth1);
        }

        // if we don't have room to taper between r1 and r2
        if ((_radius2 - _radius1) < 0.5*(_taperwidth1+_taperwidth2)) {

            // if we *really* don't have room ... taper1 by itself is too big
            // - set taper1,2 to be equal and split the r2-r1 range
            if ((_radius2 - _radius2) < 0.5*_taperwidth1) {
                _taperwidth1 = _taperwidth2 = 0.5*(_radius2 - _radius1);
                _k1 = _k2 = 1.0/(2.0*_taperwidth1);

                // if there's room for taper1, but not taper1 and 2
            } else {
                _taperwidth2 = _radius2 - _radius1 - _taperwidth1;
                _k2 = 1.0/(2.0*_taperwidth2);
            }

            _taperLo1 = _radius1 - 0.5*_taperwidth1;
            _taperHi1 = _radius1 + 0.5*_taperwidth1;
            _taperLo2 = _radius2 - 0.5*_taperwidth2;
            _taperHi2 = _radius2 + 0.5*_taperwidth2;
        }
    }

    // When called, return the throughput at the requested x,y
    // todo: replace the sinusoid taper with a band-limited
    CoordT operator() (CoordT const x, CoordT const y) const {
        CoordT const xyrad = std::sqrt(x*x + y*y);
        if ( xyrad < _taperLo1 ) {
            return 0.0;
        } else if (xyrad >= _taperLo1 && xyrad <= _taperHi1 ) {
            return 0.5*(1.0 + std::cos(  (afw::geom::TWOPI*_k1)*(xyrad - _taperHi1)));
        } else if (xyrad > _taperHi1 && xyrad <= _taperLo2 ) {
            return 1.0;
        } else if (xyrad > _taperLo2 && xyrad <= _taperHi2 ) {
            return 0.5*(1.0 + std::cos(  (afw::geom::TWOPI*_k2)*(xyrad - _taperLo2)));
        } else {
            return 0.0;
        }
    }

    CoordT getRadius1() { return _radius1; }
    CoordT getRadius2() { return _radius2; }

private:
    CoordT _radius1, _radius2;
    CoordT _taperwidth1, _taperwidth2;
    CoordT _k1, _k2;      // the angular wavenumber corresponding to a cosine with wavelength 2*taperwidth
    CoordT _taperLo1, _taperHi1;
    CoordT _taperLo2, _taperHi2;
};


template<typename CoordT>
class CircApPolar : public std::unary_function<CoordT, CoordT> {
public:
    CircApPolar(double radius, double taperwidth) : _ap(CircularAperture<CoordT>(0.0, radius, taperwidth)) {}
    CoordT operator() (double r) const {
        return r*_ap(r, 0.0);
    }
private:
    CircularAperture<CoordT> _ap;
};

/*
 * Define a Sinc functor to be integrated over for Sinc interpolation
 */
template<typename IntegrandT>
class SincAperture : public std::binary_function<IntegrandT, IntegrandT, IntegrandT> {
public:

    SincAperture(
        CircularAperture<IntegrandT> const &ap,
        int const ix,        // sinc center x
        int const iy         // sinc center y
    )
        : _ap(ap), _ix(ix), _iy(iy) {}

    IntegrandT operator() (IntegrandT const x, IntegrandT const y) const {
        double const fourierConvention = afw::geom::PI;
        double const dx = fourierConvention*(x - _ix);
        double const dy = fourierConvention*(y - _iy);
        double const fx = sinc<double>(dx);
        double const fy = sinc<double>(dy);
        return (1.0 + _ap(x, y)*fx*fy);
    }

private:
    CircularAperture<IntegrandT> const &_ap;
    double _ix, _iy;
};

class GaussPowerFunctor : public std::binary_function<double, double, double> {
public:
    GaussPowerFunctor(double sigma) : _sigma(sigma) {}

    double operator()(double kx, double ky) const {
        double const k = ::sqrt(kx*kx + ky*ky);
        double const gauss = ::exp(-0.5*k*k*_sigma*_sigma);
        return gauss*gauss;
    }
private:
    double _sigma;
};

std::pair<double, double> computeGaussLeakage(double const sigma) {

    GaussPowerFunctor gaussPower(sigma);

    double lim = afw::geom::PI;

    // total power: integrate GaussPowerFunctor -inf<x<inf, -inf<y<inf (can be done analytically)
    double powerInf = afw::geom::PI/(sigma*sigma);

    // true power: integrate GaussPowerFunctor -lim<x<lim, -lim<y<lim (must be done numerically)
    double truePower = afw::math::integrate2d(gaussPower, -lim, lim, -lim, lim, 1.0e-8);
    double trueLeak = (powerInf - truePower)/powerInf;

    // estimated power: function is circular, but coords are cartesian
    // - true power does the actual integral numerically, but we can estimate it by integrating
    //   in polar coords over lim <= radius < infinity.  The integral is analytic.
    double estLeak = ::exp(-sigma*sigma*afw::geom::PI*afw::geom::PI)/powerInf;

    return std::pair<double, double>(trueLeak, estLeak);

}

template<typename PixelT>
typename afw::image::Image<PixelT>::Ptr
calcImageRealSpace(double const rad1, double const rad2, double const taperwidth) {

    PixelT initweight = 0.0; // initialize the coeff values

    int log2   = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
    if (log2 < 3) { log2 = 3; }
    int hwid = pow(2, log2);
    int width  = 2*hwid - 1;

    int const xwidth     = width;
    int const ywidth     = width;

    int const x0 = -xwidth/2;
    int const y0 = -ywidth/2;

    // create an image to hold the coefficient image
    typename afw::image::Image<PixelT>::Ptr coeffImage =
        boost::make_shared<afw::image::Image<PixelT> >(afw::geom::ExtentI(xwidth, ywidth), initweight);
    coeffImage->setXY0(x0, y0);

    // create the aperture function object
    // determine the radius to use that makes 'radius' the effective radius of the aperture
    double tolerance = 1.0e-12;
    double dr = 1.0e-6;
    double err = 2.0*tolerance;
    double apEff = afw::geom::PI*rad2*rad2;
    double radIn = rad2;
    int maxIt = 20;
    int i = 0;
    while (err > tolerance && i < maxIt) {
        CircApPolar<double> apPolar1(radIn, taperwidth);
        CircApPolar<double> apPolar2(radIn+dr, taperwidth);
        double a1 = afw::geom::TWOPI * afw::math::integrate(apPolar1, 0.0, radIn+taperwidth, tolerance);
        double a2 = afw::geom::TWOPI * afw::math::integrate(apPolar2, 0.0, radIn+dr+taperwidth, tolerance);
        double dadr = (a2 - a1)/dr;
        double radNew = radIn - (a1 - apEff)/dadr;
        err = (a1 - apEff)/apEff;
        radIn = radNew;
        i++;
    }
    CircularAperture<double> ap(rad1, rad2, taperwidth);

    /* ******************************* */
    // integrate over the aperture

    // the limits of the integration over the sinc aperture
    double const limit = rad2 + taperwidth;
    double const x1 = -limit;
    double const x2 =  limit;
    double const y1 = -limit;
    double const y2 =  limit;

    for (int iY = y0; iY != y0 + coeffImage->getHeight(); ++iY) {
        int iX = x0;
        typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY-y0);
        for (
            typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY-y0);
            ptr != end;
            ++ptr
        ) {

            // create a sinc function in the CircularAperture at our location
            SincAperture<double> sincAp(ap, iX, iY);

            // integrate the sinc
            PixelT integral = afw::math::integrate2d(sincAp, x1, x2, y1, y2, 1.0e-8);

            // we actually integrated function+1.0 and now must subtract the excess volume
            // - just force it to zero in the corners
            double const dx = iX;
            double const dy = iY;
            *ptr = (std::sqrt(dx*dx + dy*dy) < xwidth/2) ?
                integral - (x2 - x1)*(y2 - y1) : 0.0;
            ++iX;
        }
    }


    double sum = 0.0;
    for (int iY = y0; iY != y0 + coeffImage->getHeight(); ++iY) {
        typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY-y0);
        for (
            typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY-y0);
            ptr != end; ++ptr
        ) {
            sum += *ptr;
        }
    }

#if 0                           // debugging
    coeffImage->writeFits("cimage.fits");
#endif

    return coeffImage;
}

class FftShifter {
public:
    FftShifter(int xwid) : _xwid(xwid) {}
    int shift(int x) {
        if (x >= _xwid/2) {
            return x - _xwid/2;
        } else {
            return x + _xwid/2 + 1;
        }
    }
private:
    int _xwid;
};

std::pair<double, double> rotate(double x, double y, double angle) {
    double c = ::cos(angle);
    double s = ::sin(angle);
    return std::pair<double, double>(x*c + y*s, -x*s + y*c);
}

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

template<typename PixelT>
typename afw::image::Image<PixelT>::Ptr calcImageKSpaceCplx(double const rad1, double const rad2,
                                                          double const posAng, double const ellipticity
                                                         ) {

    // we only need a half-width due to symmetry
    // make the hwid 2*rad2 so we have some buffer space and round up to the next power of 2
    int log2   = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
    if (log2 < 3) { log2 = 3; }
    int hwid = pow(2, log2);
    int wid  = 2*hwid - 1;
    int xcen = wid/2, ycen = wid/2;
    FftShifter fftshift(wid);

    boost::shared_array<std::complex<double> > cimg(new std::complex<double>[wid*wid]);
    std::complex<double> *c = cimg.get();
    // fftplan args: nx, ny, *in, *out, direction, flags
    // - done in-situ if *in == *out
    fftw_plan plan = fftw_plan_dft_2d(wid, wid,
                                      reinterpret_cast<fftw_complex*>(c),
                                      reinterpret_cast<fftw_complex*>(c),
                                      FFTW_BACKWARD, FFTW_ESTIMATE);

    // compute the k-space values and put them in the cimg array
    double const twoPiRad1 = afw::geom::TWOPI*rad1;
    double const twoPiRad2 = afw::geom::TWOPI*rad2;
    double const scale = (1.0 - ellipticity);
    for (int iY = 0; iY < wid; ++iY) {
        int const fY = fftshift.shift(iY);
        double const ky = (static_cast<double>(iY) - ycen)/wid;

        for (int iX = 0; iX < wid; ++iX) {

            int const fX = fftshift.shift(iX);
            double const kx = static_cast<double>(iX - xcen)/wid;

            // rotate
            std::pair<double, double> coo = rotate(kx, ky, posAng);
            double kxr = coo.first;
            double kyr = coo.second;
            // rescale
            double const k = ::sqrt(kxr*kxr + scale*scale*kyr*kyr);

            double const airy1 = (rad1 > 0 ? rad1*J1(twoPiRad1*k) : 0.0)/k;
            double const airy2 = rad2*J1(twoPiRad2*k)/k;
            double const airy = airy2 - airy1;

            c[fY*wid + fX] = std::complex<double>(scale*airy, 0.0);
        }
    }
    c[0] = scale*afw::geom::PI*(rad2*rad2 - rad1*rad1);

    // perform the fft and clean up after ourselves
    fftw_execute(plan);
    fftw_destroy_plan(plan);

    // put the coefficients into an image
    typename afw::image::Image<PixelT>::Ptr coeffImage =
        boost::make_shared<afw::image::Image<PixelT> >(afw::geom::ExtentI(wid, wid), 0.0);

    for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
        int iX = 0;
        typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
        for (
            typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY);\
            ptr != end;
            ++ptr
        ) {
            int fX = fftshift.shift(iX);
            int fY = fftshift.shift(iY);
            *ptr = static_cast<PixelT>(c[fY*wid + fX].real()/(wid*wid));
            iX++;
        }
    }

    // reset the origin to be the middle of the image
    coeffImage->setXY0(-wid/2, -wid/2);
    return coeffImage;
}



// I'm not sure why this doesn't work with DCT-I (REDFT00), the DCT should take advantage of symmetry
// and be much faster than the DFT.  It runs but the numbers are slightly off ...
// but I have to do it as real-to-halfcplx (R2HC) to get the correct numbers.

template<typename PixelT>
typename afw::image::Image<PixelT>::Ptr calcImageKSpaceReal(double const rad1, double const rad2) {

    // we only need a half-width due to symmertry
    // make the hwid 2*rad2 so we have some buffer space and round up to the next power of 2
    int log2   = static_cast<int>(::ceil(::log10(2.0*rad2)/log10(2.0)));
    if (log2 < 3) { log2 = 3; }
    int hwid = pow(2, log2);
    int wid  = 2*hwid - 1;
    int xcen = wid/2, ycen = wid/2;
    FftShifter fftshift(wid);

    boost::shared_array<double> cimg(new double[wid*wid]);
    double *c = cimg.get();
    // fftplan args: nx, ny, *in, *out, kindx, kindy, flags
    // - done in-situ if *in == *out
    fftw_plan plan = fftw_plan_r2r_2d(wid, wid, c, c, FFTW_R2HC, FFTW_R2HC, FFTW_ESTIMATE);

    // compute the k-space values and put them in the cimg array
    double const twoPiRad1 = afw::geom::TWOPI*rad1;
    double const twoPiRad2 = afw::geom::TWOPI*rad2;
    for (int iY = 0; iY < wid; ++iY) {

        int const fY = fftshift.shift(iY);
        double const ky = (static_cast<double>(iY) - ycen)/wid;

        for (int iX = 0; iX < wid; ++iX) {
            int const fX = fftshift.shift(iX);

            // emacs indent breaks if this isn't separte
            double const iXcen = static_cast<double>(iX - xcen);
            double const kx = iXcen/wid;

            double const k = ::sqrt(kx*kx + ky*ky);
            double const airy1 = (rad1 > 0 ? rad1*J1(twoPiRad1*k) : 0.0)/k;
            double const airy2 = rad2*J1(twoPiRad2*k)/k;
            double const airy = airy2 - airy1;
            c[fY*wid + fX] = airy;

        }
    }
    int fxy = fftshift.shift(wid/2);
    c[fxy*wid + fxy] = afw::geom::PI*(rad2*rad2 - rad1*rad1);

    // perform the fft and clean up after ourselves
    fftw_execute(plan);
    fftw_destroy_plan(plan);

    // put the coefficients into an image
    typename afw::image::Image<PixelT>::Ptr coeffImage =
        boost::make_shared<afw::image::Image<PixelT> >(afw::geom::ExtentI(wid, wid), 0.0);

    for (int iY = 0; iY != coeffImage->getHeight(); ++iY) {
        int iX = 0;
        typename afw::image::Image<PixelT>::x_iterator end = coeffImage->row_end(iY);
        for (
            typename afw::image::Image<PixelT>::x_iterator ptr = coeffImage->row_begin(iY);
            ptr != end;
            ++ptr
        ) {

            // now need to reflect the quadrant we solved to the other three
            int fX = iX < hwid ? hwid - iX - 1 : iX - hwid + 1;
            int fY = iY < hwid ? hwid - iY - 1 : iY - hwid + 1;
            *ptr = static_cast<PixelT>(c[fY*wid + fX]/(wid*wid));
            iX++;
        }
    }

    // reset the origin to be the middle of the image
    coeffImage->setXY0(-wid/2, -wid/2);
    return coeffImage;
}

} // anonymous

template<typename PixelT>
SincCoeffs<PixelT>& SincCoeffs<PixelT>::getInstance()
{
    static SincCoeffs<PixelT> instance;
    return instance;
}

template<typename PixelT>
void SincCoeffs<PixelT>::cache(float r1, float r2)
{
    if (r1 < 0.0 || r2 < r1) {
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                          (boost::format("Invalid r1,r2 = %f,%f") % r1 % r2).str());
    }
    double const innerFactor = r1/r2;
    afw::geom::ellipses::Axes axes(r2, r2, 0.0);
    if (!getInstance()._lookup(axes, innerFactor)) {
        PTR(typename SincCoeffs<PixelT>::CoeffT) coeff = calculate(axes, innerFactor);
        coeff->markPersistent();
        getInstance()._cache[r2][innerFactor] = coeff;
    }
}

template<typename PixelT>
CONST_PTR(typename SincCoeffs<PixelT>::CoeffT)
SincCoeffs<PixelT>::get(afw::geom::ellipses::Axes const& axes, float const innerFactor)
{
    CONST_PTR(CoeffT) coeff = getInstance()._lookup(axes, innerFactor);
    return coeff ? coeff : calculate(axes, innerFactor);
}

template<typename PixelT>
CONST_PTR(typename SincCoeffs<PixelT>::CoeffT)
SincCoeffs<PixelT>::_lookup(afw::geom::ellipses::Axes const& axes, double const innerFactor) const
{
    if (innerFactor < 0.0 || innerFactor > 1.0) {
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                          (boost::format("innerFactor = %f is not between 0 and 1") % innerFactor).str());
    }

    CONST_PTR(typename SincCoeffs<PixelT>::CoeffT) const null = CONST_PTR(SincCoeffs<PixelT>::CoeffT)();

    // We only cache circular apertures
    if (!FuzzyCompare<float>().isEqual(axes.getA(), axes.getB())) {
        return null;
    }
    typename CoeffMapMap::const_iterator iter1 = _cache.find(axes.getA());
    if (iter1 == _cache.end()) {
        return null;
    }
    typename CoeffMap::const_iterator iter2 = iter1->second.find(innerFactor);
    return (iter2 == iter1->second.end()) ? null : iter2->second;
}

template<typename PixelT>
PTR(typename SincCoeffs<PixelT>::CoeffT)
SincCoeffs<PixelT>::calculate(afw::geom::ellipses::Axes const& axes, double const innerFactor)
{
    if (innerFactor < 0.0 || innerFactor > 1.0) {
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                          (boost::format("innerFactor = %f is not between 0 and 1") % innerFactor).str());
    }

    // Kspace-real is fastest, but only slightly faster than kspace cplx
    // but real won't work for elliptical apertures due to symmetries assumed for real transform

    double const rad1 = axes.getA() * innerFactor;
    double const rad2 = axes.getA();
    // if there's no angle and no ellipticity
    if (FuzzyCompare<float>().isEqual(axes.getA(), axes.getB())) {
        // here we call the real transform
        return calcImageKSpaceReal<PixelT>(rad1, rad2);
    } else {
        // here we call the complex transform
        double const ellipticity = 1.0 - axes.getB()/axes.getA();
        return calcImageKSpaceCplx<PixelT>(rad1, rad2, axes.getTheta(), ellipticity);
    }
}

template class SincCoeffs<float>;
template class SincCoeffs<double>;

}}} // namespace lsst::meas::base