1# This file is part of the Astrometry.net suite.
2# Licensed under a 3-clause BSD style license - see LICENSE
3from __future__ import print_function
4import matplotlib
5matplotlib.use('Agg')
6
7import sys
8from optparse import *
9
10import numpy as np
11from pylab import *
12from numpy import *
13#from astrometry.util.sip import *
14from astrometry.util.util import *
15
16def plot_distortions(wcsfn, ex=1, ngridx=10, ngridy=10, stepx=10, stepy=10):
17    wcs = Sip(wcsfn)
18    W,H = wcs.wcstan.imagew, wcs.wcstan.imageh
19
20    xgrid = np.linspace(0, W, ngridx)
21    ygrid = np.linspace(0, H, ngridy)
22    X = np.linspace(0, W, int(ceil(W/stepx)))
23    Y = np.linspace(0, H, int(ceil(H/stepy)))
24
25    xlo,xhi,ylo,yhi = 0,W,0,H
26
27    for x in xgrid:
28        DX,DY = [],[]
29        xx,yy = [],[]
30        UX,UY = [],[]
31        for y in Y:
32            dx,dy = wcs.get_distortion(x, y)
33            ux,uy = wcs.get_undistortion(dx, dy)
34            print('x,y', (x,y), 'dx,dy', (dx,dy), 'ux,uy', (ux,uy))
35            xx.append(x)
36            yy.append(y)
37            DX.append(dx)
38            DY.append(dy)
39            UX.append(ux)
40            UY.append(uy)
41        DX = np.array(DX)
42        DY = np.array(DY)
43        UX = np.array(UX)
44        UY = np.array(UY)
45        xx = np.array(xx)
46        yy = np.array(yy)
47        EX = DX + ex * (DX - xx)
48        EY = DY + ex * (DY - yy)
49        plot(xx, yy, 'k-', alpha=0.5)
50        plot(EX, EY, 'r-')
51        plot(UX, UY, 'b-', alpha=0.5)
52        xlo = min(xlo, min(EX))
53        xhi = max(xhi, max(EX))
54        ylo = min(ylo, min(EY))
55        yhi = max(yhi, max(EY))
56
57    for y in ygrid:
58        DX,DY = [],[]
59        xx,yy = [],[]
60        UX,UY = [],[]
61        for x in X:
62            dx,dy = wcs.get_distortion(x, y)
63            ux,uy = wcs.get_undistortion(dx, dy)
64            DX.append(dx)
65            DY.append(dy)
66            UX.append(ux)
67            UY.append(uy)
68            xx.append(x)
69            yy.append(y)
70        DX = np.array(DX)
71        DY = np.array(DY)
72        xx = np.array(xx)
73        yy = np.array(yy)
74        EX = DX + ex * (DX - xx)
75        EY = DY + ex * (DY - yy)
76        plot(xx, yy, 'k-', alpha=0.5)
77        plot(EX, EY, 'r-')
78        plot(UX, UY, 'b-', alpha=0.5)
79        xlo = min(xlo, min(EX))
80        xhi = max(xhi, max(EX))
81        ylo = min(ylo, min(EY))
82        yhi = max(yhi, max(EY))
83
84    plot([wcs.wcstan.crpix[0]], [wcs.wcstan.crpix[1]], 'rx')
85
86    #axis([0, W, 0, H])
87    axis('scaled')
88    axis([xlo,xhi,ylo,yhi])
89    #axis('tight')
90
91if __name__ == '__main__':
92    parser = OptionParser(usage='%prog [options] <wcs-filename> <plot-filename>')
93    parser.add_option('-e', '--ex', '--exaggerate', dest='ex', type='float', help='Exaggerate the distortion by this factor')
94    #parser.add_option('-s', '--scale', dest='scale', type='float', help='Scale the
95    parser.add_option('-n', dest='nsteps', type='int', help='Number of grid lines to plot')
96
97    parser.set_defaults(ex=1.)
98
99    opt,args = parser.parse_args()
100
101    if len(args) != 2:
102        parser.print_help()
103        sys.exit(-1)
104
105    wcsfn = args[0]
106    outfn = args[1]
107
108    args = {}
109    if opt.ex is not None:
110        args['ex'] = opt.ex
111    if opt.nsteps is not None:
112        args['ngridx'] = opt.nsteps
113        args['ngridy'] = opt.nsteps
114
115    clf()
116    plot_distortions(wcsfn, **args)
117    tt = 'SIP distortions: %s' % wcsfn
118    if opt.ex != 1:
119        tt += ' (exaggerated by %g)' % opt.ex
120    title(tt)
121    savefig(outfn)
122
123