How to use the tractor.ducks.ImageCalibration function in tractor

from __future__ import print_function
import numpy as np
import scipy.interpolate as interp
from tractor.utils import ParamList
from tractor import ducks


class SplineSky(ParamList, ducks.ImageCalibration):

    @staticmethod
    def BlantonMethod(image, mask, gridsize, estimator=np.median):
        '''
        mask: True to use pixel.  None for no masking.
        '''
        H, W = image.shape
        halfbox = gridsize // 2

        nx = int(np.round(W // float(halfbox))) + 2
        x0 = int((W - (nx - 2) * halfbox) // 2)
        xgrid = x0 + (halfbox * (np.arange(nx) - 0.5)).astype(int)

        ny = int(np.round(H // float(halfbox))) + 2
        y0 = int((H - (ny - 2) * halfbox) // 2)
        ygrid = y0 + (halfbox * (np.arange(ny) - 0.5)).astype(int)

from tractor.utils import BaseParams, ScalarParam
from tractor.patch import Patch
from tractor import ducks

class NullSky(BaseParams, ducks.Sky):
    '''
    A Sky implementation that does nothing; the background level is
    zero.
    '''
    pass
    
class ConstantSky(ScalarParam, ducks.ImageCalibration):
    '''
    A simple constant sky level across the whole image.

    This sky object has one parameter, the constant level.
    
    The sky level is specified in the same units as the image
    ("counts").
    '''
    def getParamDerivatives(self, tractor, img, srcs):
        import numpy as np
        p = Patch(0, 0, np.ones_like(img.getImage()))
        p.setName('dsky')
        return [p]
    def addTo(self, img, scale=1.):
        if self.val == 0:
            return

# print 'Fitting:'
        # tr.printThawedParams()
        tim.modelMinval = approx
        alphas = [0.1, 0.3, 1.0]
        for step in range(50):
            dlnp, X, alpha = tr.optimize(shared_params=False, alphas=alphas,
                                         damp=damp)
            # print 'dlnp', dlnp, 'alpha', alpha
            # print 'X', X
            if dlnp < 1e-6:
                break
            # print 'psf', psf
        return psf


class NCircularGaussianPSF(MultiParams, ducks.ImageCalibration):
    '''
    A PSF model using N concentric, circular Gaussians.
    '''

    def __init__(self, sigmas, weights):
        '''
        Creates a new N-Gaussian (concentric, isotropic) PSF.

        sigmas: (list of floats) standard deviations of the components

        weights: (list of floats) relative weights of the components;
        given two components with weights 0.9 and 0.1, the total mass
        due to the second component will be 0.1.  These values will be
        normalized so that the total mass of the PSF is 1.0.

        eg,   NCircularGaussianPSF([1.5, 4.0], [0.8, 0.2])

#         for ii in range(len(axes)-1):
#             a = m.ifft(a, s[ii], axes[ii], overwrite_x=ovr_x)
#             ovr_x = True
#     a = m.irfft_numpy(a, n = s[-1], axis=la)
#     return a
# m.irfftn_numpy  = irfftn_numpy

from tractor.utils import MultiParams, getClassName
from tractor.psf import GaussianMixturePSF, PixelizedPSF
from tractor import mixture_profiles as mp
from tractor import ducks

from astrometry.util.fits import fits_table


class VaryingGaussianPSF(MultiParams, ducks.ImageCalibration):
    '''
    A mixture-of-Gaussians (MoG) PSF with spatial variation,
    represented as a spline(x,y) in each of the MoG parameters.

    This is a base class -- subclassers must implement "instantiateAt"
    '''

    def __init__(self, W, H, nx=11, ny=11, K=3, psfClass=GaussianMixturePSF):
        '''
        W,H: image size (for the image where this PSF lives)
        nx,ny: number of sample points for the spline fit
        K: number of components in the Gaussian mixture
        '''
        # HACK -- no args?
        super(VaryingGaussianPSF, self).__init__()
        self.psfclass = psfClass

# print 'PixelizedPSF FFT size', sz
        if sz in self.fftcache:
            return self.fftcache[sz]

        pad, cx, cy = self._padInImage(sz, sz)
        # cx,cy: coordinate of the PSF center in *pad*
        P = np.fft.rfft2(pad)
        P = P.astype(np.complex64)
        pH, pW = pad.shape
        v = np.fft.rfftfreq(pW)
        w = np.fft.fftfreq(pH)
        rtn = P, (cx, cy), (pH, pW), (v, w)
        self.fftcache[sz] = rtn
        return rtn

class GaussianMixturePSF(MogParams, ducks.ImageCalibration):
    '''
    A PSF model that is a mixture of general 2-D Gaussians
    (characterized by amplitude, mean, covariance)
    '''

    def __init__(self, *args):
        '''
        GaussianMixturePSF(amp, mean, var)

        or

        GaussianMixturePSF(a0,a1,a2, mx0,my0,mx1,my1,mx2,my2,
                           vxx0,vyy0,vxy0, vxx1,vyy1,vxy1, vxx2,vyy2,vxy2)

        amp:  np array (size K) of Gaussian amplitudes
        mean: np array (size K,2) of means

k = prefix + 'A%i' % i
            if not k in hdr:
                break
            args.append(hdr.get(k))
        obj = clazz(*args)
        params = []
        for i in range(100):
            k = prefix + 'P%i' % i
            if not k in hdr:
                break
            params.append(hdr.get(k))
        obj.setAllParams(params)
        return obj


class Sky(ImageCalibration, Params):
    '''
    Duck-type definition for a sky model.
    '''

    def getParamDerivatives(self, tractor, img, srcs):
        '''
        Returns [ Patch, Patch, ... ], of length numberOfParams(),
        containing the derivatives in the given `Image` for each
        parameter.
        '''
        return []

    def addTo(self, mod, scale=1.):
        '''
        Add the sky to the input synthetic image `mod`, a 2-D numpy
        array.

'''
        Returns the local pixel scale at the given *x*,*y* pixel coords,
        in *arcseconds* per pixel.
        '''
        import numpy as np
        return 3600. * np.sqrt(np.abs(np.linalg.det(self.cdAtPixel(x, y))))

    def shifted(self, dx, dy):
        '''
        Returns a new WCS object appropriate for the subimage starting
        at (dx,dy) with respect to the current WCS origin.
        '''
        return None


class PSF(ImageCalibration, Params):
    '''
    Duck-type definition of a point-spread function.
    '''

    def getPointSourcePatch(self, px, py, minval=0., modelMask=None):
        '''
        Returns a `Patch`, a rendering of a point source at the given
        pixel coordinates.

        The returned `Patch` should have unit "counts".

        *minval* says that we are willing to accept an approximation
        such that pixels with counts < minval can be omitted.

        *modelMask* describes the pixels to be evaluated.  If the
         *modelMask* includes a pixel-by-pixel mask, this overrides

(Returns the constant ``CD`` matrix elements)
        '''
        ok, ra0, dec0 = self.wcs.pixelxy2radec(
            x + 1. + self.x0, y + 1. + self.y0)
        ok, ra1, dec1 = self.wcs.pixelxy2radec(
            x + 2. + self.x0, y + 1. + self.y0)
        ok, ra2, dec2 = self.wcs.pixelxy2radec(
            x + 1. + self.x0, y + 2. + self.y0)

        cosdec = np.cos(np.deg2rad(dec0))

        return np.array([[(ra1 - ra0) * cosdec, (ra2 - ra0) * cosdec],
                         [dec1 - dec0,        dec2 - dec0]])


class ConstantFitsWcs(ParamList, ducks.ImageCalibration):
    '''
    A WCS implementation that wraps a FITS WCS object (with a pixel
    offset).
    '''

    def __init__(self, wcs):
        '''
        Creates a new ``ConstantFitsWcs`` given an underlying WCS object.
        '''
        self.x0 = 0
        self.y0 = 0
        super(ConstantFitsWcs, self).__init__()
        self.wcs = wcs

        cd = self.wcs.get_cd()
        self.cd = np.array([[cd[0], cd[1]], [cd[2], cd[3]]])

Returns the local pixel scale at the given *x*,*y* pixel coords,
        in *arcseconds* per pixel.
        '''
        return self.pixscale

    def toStandardFitsHeader(self, hdr):
        if self.x0 != 0 or self.y0 != 0:
            wcs = self.wcs.get_subimage(self.x0, self.y0,
                                        self.wcs.get_width() - self.x0,
                                        self.wcs.get_height() - self.y0)
        else:
            wcs = self.wcs
        wcs.add_to_header(hdr)


class TanWcs(ConstantFitsWcs, ducks.ImageCalibration):
    '''
    The WCS object must be an astrometry.util.util.Tan object, or a
    convincingly quacking duck.

    You can also give it a filename and HDU.
    '''

    @staticmethod
    def getNamedParams():
        return dict(crval1=0, crval2=1, crpix1=2, crpix2=3,
                    cd1_1=4, cd1_2=5, cd2_1=6, cd2_2=7)

    def __init__(self, wcs, hdu=0):
        '''
        Creates a new ``FitsWcs`` given a :class:`~astrometry.util.util.Tan`
        object.  To create one of these from a filename and FITS HDU extension,

def __init__(self, band):
        self.band = band
        BaseParams.__init__(self)

    def copy(self):
        return FluxesPhotoCal(self.band)

    def brightnessToCounts(self, brightness):
        flux = brightness.getFlux(self.band)
        return flux

    def __str__(self):
        return 'FluxesPhotoCal(band=%s)' % (self.band)


class MagsPhotoCal(ParamList, ducks.ImageCalibration):
    '''
    A `PhotoCal` implementation to be used with zeropoint-calibrated
    `Mags`.
    '''

    def __init__(self, band, zeropoint):
        '''
        Create a new ``MagsPhotoCal`` object with a *zeropoint* in a
        *band*.

        The ``Mags`` objects you use must have *band* as one of their
        available bands.
        '''
        self.band = band
        # MAGIC
        self.maxmag = 50.