1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 #include <math.h>
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <string.h>
10 #include <assert.h>
11
12 #include "tpv.h"
13 #include "starutil.h"
14 #include "mathutil.h"
15
16 /*******************************
17
18 NOTE, this does not work yet!!!
19
20 I copied this from sip.c, which applies the distortion in PIXEL space,
21 while TPV applies the distortion in IWC space.
22
23
24
25 The TPV projection is evaluated as follows.
26
27 Compute the first order standard coordinates xi and eta from the
28 linear part of the solution stored in CRPIX and the CD matrix.
29
30 xi = CD1_1 * (x - CRPIX1) + CD1_2 * (y - CRPIX2)
31 eta = CD2_1 * (x - CRPIX1) + CD2_2 * (y - CRPIX2)
32
33 Apply the distortion transformation using the coefficients in the
34 PV keywords as described below.
35
36 xi' = f_xi (xi, eta)
37 eta' = f_eta (xi, eta)
38
39 Apply the tangent plane projection to xi' and eta' as described in
40 Calabretta and Greisen . The reference tangent point given by the
41 CRVAL values lead to the final RA and DEC in degrees. Note that
42 the units of xi, eta, f_xi, and f_eta are also degrees.
43
44 ********************************/
45
46
tpv_xyz2pixelxy(const tpv_t * tpv,double x,double y,double z,double * px,double * py)47 anbool tpv_xyz2pixelxy(const tpv_t* tpv, double x, double y, double z, double *px, double *py) {
48 double xyz[3];
49 xyz[0] = x;
50 xyz[1] = y;
51 xyz[2] = z;
52 return tpv_xyzarr2pixelxy(tpv, xyz, px, py);
53 }
54
tpv_wrap_tan(const tan_t * tan,tpv_t * tpv)55 void tpv_wrap_tan(const tan_t* tan, tpv_t* tpv) {
56 memset(tpv, 0, sizeof(tpv_t));
57 memcpy(&(tpv->wcstan), tan, sizeof(tan_t));
58 }
59
tpv_get_crval(const tpv_t * tpv,double * ra,double * dec)60 void tpv_get_crval(const tpv_t* tpv, double* ra, double* dec) {
61 *ra = tpv->wcstan.crval[0];
62 *dec = tpv->wcstan.crval[1];
63 }
64
tpv_imagew(tpv_t * tpv)65 double tpv_imagew(tpv_t* tpv) {
66 assert(tpv);
67 return tpv->wcstan.imagew;
68 }
tpv_imageh(tpv_t * tpv)69 double tpv_imageh(tpv_t* tpv) {
70 assert(tpv);
71 return tpv->wcstan.imageh;
72 }
tpv_create()73 tpv_t* tpv_create() {
74 tpv_t* tpv = calloc(1, sizeof(tpv_t));
75
76 tpv->wcstan.cd[0][0] = 1;
77 tpv->wcstan.cd[0][1] = 0;
78 tpv->wcstan.cd[1][0] = 0;
79 tpv->wcstan.cd[1][1] = 1;
80
81 return tpv;
82 }
83
tpv_free(tpv_t * tpv)84 void tpv_free(tpv_t* tpv) {
85 free(tpv);
86 }
87
tpv_copy(tpv_t * dest,const tpv_t * src)88 void tpv_copy(tpv_t* dest, const tpv_t* src) {
89 memcpy(dest, src, sizeof(tpv_t));
90 }
91
tpv_distortion(const tpv_t * tpv,double x,double y,double * X,double * Y)92 static void tpv_distortion(const tpv_t* tpv, double x, double y,
93 double* X, double* Y) {
94 // Get pixel coordinates relative to reference pixel
95 double u = x - tpv->wcstan.crpix[0];
96 double v = y - tpv->wcstan.crpix[1];
97 tpv_calc_distortion(tpv, u, v, X, Y);
98 *X += tpv->wcstan.crpix[0];
99 *Y += tpv->wcstan.crpix[1];
100 }
101
102 // Pixels to RA,Dec in degrees.
tpv_pixelxy2radec(const tpv_t * tpv,double px,double py,double * ra,double * dec)103 void tpv_pixelxy2radec(const tpv_t* tpv, double px, double py,
104 double *ra, double *dec) {
105 double U, V;
106 tpv_distortion(tpv, px, py, &U, &V);
107 // Run a normal TAN conversion on the distorted pixel coords.
108 tan_pixelxy2radec(&(tpv->wcstan), U, V, ra, dec);
109 }
110
111 // Pixels to Intermediate World Coordinates in degrees.
tpv_pixelxy2iwc(const tpv_t * tpv,double px,double py,double * iwcx,double * iwcy)112 void tpv_pixelxy2iwc(const tpv_t* tpv, double px, double py,
113 double *iwcx, double* iwcy) {
114 double U, V;
115 tpv_distortion(tpv, px, py, &U, &V);
116 // Run a normal TAN conversion on the distorted pixel coords.
117 tan_pixelxy2iwc(&(tpv->wcstan), U, V, iwcx, iwcy);
118 }
119
120 // Pixels to XYZ unit vector.
tpv_pixelxy2xyzarr(const tpv_t * tpv,double px,double py,double * xyz)121 void tpv_pixelxy2xyzarr(const tpv_t* tpv, double px, double py, double *xyz) {
122 double U, V;
123 tpv_distortion(tpv, px, py, &U, &V);
124 // Run a normal TAN conversion on the distorted pixel coords.
125 tan_pixelxy2xyzarr(&(tpv->wcstan), U, V, xyz);
126 }
127
tpv_iwc2pixelxy(const tpv_t * tpv,double u,double v,double * px,double * py)128 void tpv_iwc2pixelxy(const tpv_t* tpv, double u, double v,
129 double *px, double* py) {
130 double x,y;
131 tan_iwc2pixelxy(&(tpv->wcstan), u, v, &x, &y);
132 tpv_pixel_undistortion(tpv, x, y, px, py);
133 }
134
135 // RA,Dec in degrees to Pixels.
tpv_radec2pixelxy(const tpv_t * tpv,double ra,double dec,double * px,double * py)136 anbool tpv_radec2pixelxy(const tpv_t* tpv, double ra, double dec, double *px, double *py) {
137 double x,y;
138 if (!tan_radec2pixelxy(&(tpv->wcstan), ra, dec, &x, &y))
139 return FALSE;
140 tpv_pixel_undistortion(tpv, x, y, px, py);
141 return TRUE;
142 }
143
tpv_iwc2radec(const tpv_t * tpv,double x,double y,double * p_ra,double * p_dec)144 void tpv_iwc2radec(const tpv_t* tpv, double x, double y, double *p_ra, double *p_dec) {
145 tan_iwc2radec(&(tpv->wcstan), x, y, p_ra, p_dec);
146 }
147
tpv_xyzarr2pixelxy(const tpv_t * tpv,const double * xyz,double * px,double * py)148 anbool tpv_xyzarr2pixelxy(const tpv_t* tpv, const double* xyz, double *px, double *py) {
149 double ra, dec;
150 xyzarr2radecdeg(xyz, &ra, &dec);
151 return tpv_radec2pixelxy(tpv, ra, dec, px, py);
152 }
153
154
tpv_xyzarr2iwc(const tpv_t * tpv,const double * xyz,double * iwcx,double * iwcy)155 anbool tpv_xyzarr2iwc(const tpv_t* tpv, const double* xyz,
156 double* iwcx, double* iwcy) {
157 return tan_xyzarr2iwc(&(tpv->wcstan), xyz, iwcx, iwcy);
158 }
tpv_radec2iwc(const tpv_t * tpv,double ra,double dec,double * iwcx,double * iwcy)159 anbool tpv_radec2iwc(const tpv_t* tpv, double ra, double dec,
160 double* iwcx, double* iwcy) {
161 return tan_radec2iwc(&(tpv->wcstan), ra, dec, iwcx, iwcy);
162 }
163
tpv_calc_distortion(const tpv_t * tpv,double xi,double eta,double * XI,double * ETA)164 void tpv_calc_distortion(const tpv_t* tpv, double xi, double eta,
165 double* XI, double *ETA) {
166 // Do TPV distortion (in intermediate world coords)
167 // xi,eta -> xi',eta'
168
169 /**
170 From http://iraf.noao.edu/projects/mosaic/tpv.html
171
172 p = PV1_
173
174 xi' = p0 +
175 p1 * xi + p2 * eta + p3 * r +
176 p4 * xi^2 + p5 * xi * eta + p6 * eta^2 +
177 p7 * xi^3 + p8 * xi^2 * eta + p9 * xi * eta^2 +
178 p10 * eta^3 + p11 * r^3 +
179 p12 * xi^4 + p13 * xi^3 * eta + p14 * xi^2 * eta^2 +
180 p15 * xi * eta^3 + p16 * eta^4 +
181 p17 * xi^5 + p18 * xi^4 * eta + p19 * xi^3 * eta^2 +
182 p20 * xi^2 * eta^3 + p21 * xi * eta^4 + p22 * eta^5 + p23 * r^5 +
183 p24 * xi^6 + p25 * xi^5 * eta + p26 * xi^4 * eta^2 +
184 p27 * xi^3 * eta^3 + p28 * xi^2 * eta^4 + p29 * xi * eta^5 +
185 p30 * eta^6
186 p31 * xi^7 + p32 * xi^6 * eta + p33 * xi^5 * eta^2 +
187 p34 * xi^4 * eta^3 + p35 * xi^3 * eta^4 + p36 * xi^2 * eta^5 +
188 p37 * xi * eta^6 + p38 * eta^7 + p39 * r^7
189
190 p = PV2_
191 eta' = p0 +
192 p1 * eta + p2 * xi + p3 * r +
193 p4 * eta^2 + p5 * eta * xi + p6 * xi^2 +
194 p7 * eta^3 + p8 * eta^2 * xi + p9 * eta * xi^2 + p10 * xi^3 +
195 p11 * r^3 +
196 p12 * eta^4 + p13 * eta^3 * xi + p14 * eta^2 * xi^2 +
197 p15 * eta * xi^3 + p16 * xi^4 +
198 p17 * eta^5 + p18 * eta^4 * xi + p19 * eta^3 * xi^2 +
199 p20 * eta^2 * xi^3 + p21 * eta * xi^4 + p22 * xi^5 +
200 p23 * r^5 +
201 p24 * eta^6 + p25 * eta^5 * xi + p26 * eta^4 * xi^2 +
202 p27 * eta^3 * xi^3 + p28 * eta^2 * xi^4 + p29 * eta * xi^5 +
203 p30 * xi^6
204 p31 * eta^7 + p32 * eta^6 * xi + p33 * eta^5 * xi^2 +
205 p34 * eta^4 * xi^3 + p35 * eta^3 * xi^4 + p36 * eta^2 * xi^5 +
206 p37 * eta * xi^6 + p38 * xi^7 + p39 * r^7
207
208 Note the "cross-over" -- the xi' powers are in terms of xi,eta
209 while the eta' powers are in terms of eta,xi.
210 */
211
212 // 1 x y r x2 xy y2 x3 x2y xy2 y3 r3 x4 x3y x2y2 xy3 y4
213 // x5 x4y x3y2 x2y3 xy4 y5 r5 x6 x5y x4y2, x3y3 x2y4 xy5 y6
214 // x7 x6y x5y2 x4y3 x3y4 x2y5 xy6 y7 r7
215 int xp[] = {
216 0,
217 1, 0, 0,
218 2, 1, 0,
219 3, 2, 1, 0, 0,
220 4, 3, 2, 1, 0,
221 5, 4, 3, 2, 1, 0, 0,
222 6, 5, 4, 3, 2, 1, 0,
223 7, 6, 5, 4, 3, 2, 1, 0, 0};
224 int yp[] = {
225 0,
226 0, 1, 0,
227 0, 1, 2,
228 0, 1, 2, 3, 0,
229 0, 1, 2, 3, 4,
230 0, 1, 2, 3, 4, 5, 0,
231 0, 1, 2, 3, 4, 5, 6,
232 0, 1, 2, 3, 4, 5, 6, 7, 0};
233 int rp[] = {
234 0,
235 0, 0, 1,
236 0, 0, 0,
237 0, 0, 0, 0, 3,
238 0, 0, 0, 0, 0,
239 0, 0, 0, 0, 0, 0, 5,
240 0, 0, 0, 0, 0, 0, 0,
241 0, 0, 0, 0, 0, 0, 0, 0, 7};
242 double xipows[8];
243 double etapows[8];
244 double rpows[8];
245 double r;
246 int i;
247 double px, py;
248
249 r = sqrt(xi*xi + eta*eta);
250 xipows[0] = etapows[0] = rpows[0] = 1.0;
251 for (i=1; i<sizeof(xipows)/sizeof(double); i++) {
252 xipows[i] = xipows[i-1]*xi;
253 etapows[i] = etapows[i-1]*eta;
254 rpows[i] = rpows[i-1]*r;
255 }
256 px = py = 0;
257 for (i=0; i<sizeof(xp)/sizeof(int); i++) {
258 px += tpv->pv1[i] * xipows[xp[i]] * etapows[yp[i]]
259 * rpows[rp[i]];
260 // here's the "cross-over" mentioned above
261 py += tpv->pv2[i] * etapows[xp[i]] * xipows[yp[i]]
262 * rpows[rp[i]];
263 }
264 *XI = px;
265 *ETA = py;
266 }
267
tpv_pixel_distortion(const tpv_t * tpv,double x,double y,double * X,double * Y)268 void tpv_pixel_distortion(const tpv_t* tpv, double x, double y, double* X, double *Y) {
269 tpv_distortion(tpv, x, y, X, Y);
270 }
271
tpv_pixel_undistortion(const tpv_t * tpv,double x,double y,double * X,double * Y)272 void tpv_pixel_undistortion(const tpv_t* tpv, double x, double y, double* X, double *Y) {
273 // Get pixel coordinates relative to reference pixel
274 double u = x - tpv->wcstan.crpix[0];
275 double v = y - tpv->wcstan.crpix[1];
276 tpv_calc_inv_distortion(tpv, u, v, X, Y);
277 *X += tpv->wcstan.crpix[0];
278 *Y += tpv->wcstan.crpix[1];
279 }
280
tpv_calc_inv_distortion(const tpv_t * tpv,double U,double V,double * u,double * v)281 void tpv_calc_inv_distortion(const tpv_t* tpv, double U, double V, double* u, double *v)
282 {
283 assert(0);
284 printf("tpv_calc_inv_distortion not implemented!\n");
285 exit(-1);
286 }
287
tpv_det_cd(const tpv_t * tpv)288 double tpv_det_cd(const tpv_t* tpv) {
289 return tan_det_cd(&(tpv->wcstan));
290 }
291
292 // returns pixel scale in arcseconds (NOT arcsec^2)
tpv_pixel_scale(const tpv_t * tpv)293 double tpv_pixel_scale(const tpv_t* tpv) {
294 return tan_pixel_scale(&(tpv->wcstan));
295 }
296
tpv_print_to(const tpv_t * tpv,FILE * f)297 void tpv_print_to(const tpv_t* tpv, FILE* f) {
298 double det,pixsc;
299
300 print_to(&(tpv->wcstan), f, "TPV");
301
302 /*
303 if (tpv->a_order > 0) {
304 int p, q;
305 for (p=0; p<=tpv->a_order; p++) {
306 fprintf(f, (p ? " " : " A = "));
307 for (q=0; q<=tpv->a_order; q++)
308 if (p+q <= tpv->a_order)
309 //fprintf(f,"a%d%d=%le\n", p,q,tpv->a[p][q]);
310 fprintf(f,"%12.5g", tpv->a[p][q]);
311 fprintf(f,"\n");
312 }
313 }
314 if (tpv->b_order > 0) {
315 int p, q;
316 for (p=0; p<=tpv->b_order; p++) {
317 fprintf(f, (p ? " " : " B = "));
318 for (q=0; q<=tpv->b_order; q++)
319 if (p+q <= tpv->a_order)
320 fprintf(f,"%12.5g", tpv->b[p][q]);
321 //if (p+q <= tpv->b_order && p+q > 0)
322 //fprintf(f,"b%d%d=%le\n", p,q,tpv->b[p][q]);
323 fprintf(f,"\n");
324 }
325 }
326 det = tpv_det_cd(tpv);
327 pixsc = 3600*sqrt(fabs(det));
328 //fprintf(f," det(CD)=%g\n", det);
329 fprintf(f," sqrt(det(CD))=%g [arcsec]\n", pixsc);
330 //fprintf(f,"\n");
331 */
332 }
333
tpv_print(const tpv_t * tpv)334 void tpv_print(const tpv_t* tpv) {
335 tpv_print_to(tpv, stderr);
336 }
337
tpv_get_orientation(const tpv_t * tpv)338 double tpv_get_orientation(const tpv_t* tpv) {
339 return tan_get_orientation(&(tpv->wcstan));
340 }
341