1 /*********************************************************************
2 wcsdistortion -- Manipulation of distortions in WCS.
3 This is part of GNU Astronomy Utilities (Gnuastro) package.
4
5 Original author:
6 Sachin Kumar Singh <sachinkumarsingh092@gmail.com>
7 Contributing author(s):
8 Mohammad Akhlaghi <mohammad@akhlaghi.org>
9 Copyright (C) 2020-2021, Free Software Foundation, Inc.
10
11 Gnuastro is free software: you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation, either version 3 of the License, or (at your
14 option) any later version.
15
16 Gnuastro is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 General Public License for more details.
20
21 You should have received a copy of the GNU General Public License
22 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
23 **********************************************************************/
24 #include <config.h>
25
26 #include <errno.h>
27 #include <error.h>
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <string.h>
31 #include <unistd.h>
32
33 #include <wcslib/wcslib.h>
34
35 #include <gsl/gsl_linalg.h>
36 #include <gsl/gsl_vector.h>
37 #include <gsl/gsl_matrix.h>
38 #include <gsl/gsl_multifit.h>
39
40 #include <gnuastro/wcs.h>
41 #include <gnuastro/fits.h>
42
43 #include <gnuastro-internal/wcsdistortion.h>
44
45
46
47
48
49 /* Internally used macro(s) to help in the processing */
50 #define wcsdistortion_max(a,b) \
51 ({ __typeof__ (a) _a = (a); \
52 __typeof__ (b) _b = (b); \
53 _a > _b ? _a : _b; })
54
55 /* Declarations to avoid ordering problems. */
56 static void
57 wcsdistortion_calc_tpveq(struct wcsprm *wcs, double cd[2][2],
58 double tpvu[8][8], double tpvv[8][8]);
59 static double
60 wcsdistortion_calcsip(size_t axis, size_t m, size_t n, double tpvu[8][8],
61 double tpvv[8][8]);
62
63
64
65
66
67
68
69 /**************************************************************/
70 /********** Reading utilities ************/
71 /**************************************************************/
72 /* Extract the required pv parameters from wcs. */
73 static void
wcsdistortion_get_tpvparams(struct wcsprm * wcs,double cd[2][2],double * pv1,double * pv2)74 wcsdistortion_get_tpvparams(struct wcsprm *wcs, double cd[2][2],
75 double *pv1, double *pv2)
76 {
77 const char *cp;
78 double *temp_cd=NULL;
79 struct dpkey *keyp=NULL;
80 size_t i, j, k=0, index=0;
81 struct disprm *disseq=NULL;
82
83 /* Make sure a WCS structure is actually given. */
84 if(wcs==NULL)
85 error(EXIT_FAILURE, 0, "%s: input WCS structure is NULL", __func__);
86
87 /* For easy reading. */
88 disseq=wcs->lin.disseq;
89 keyp=disseq->dp;
90
91 /* Fill the 2-times allocated CD array (cd[][]). Note that the required
92 CD matix is extracted using the `gal_wcs_wrap_matrix` as a single
93 allocated array (temp_cd[]), that is then used to fill cd[][]. */
94 temp_cd=gal_wcs_warp_matrix(wcs);
95 for(i=0; i<2; ++i)
96 for(j=0; j<2; ++j)
97 {
98 /* If a value is present store it, else store 0.*/
99 if(temp_cd[k] != 0) cd[i][j]=temp_cd[k++];
100 else cd[i][j]=0;
101
102 /* For a check:
103 printf("CD%ld_%ld\t%.10E\n", i+1, j+1, cd[i][j]);
104 */
105 }
106
107 for(i=0; i<disseq->ndp; ++i, ++keyp)
108 {
109 /* For axis1. */
110 if (keyp->j == 1)
111 {
112 if ( keyp->field[0] == '\0' ) continue;
113 cp = strchr(keyp->field, '.') + 1;
114 if (strncmp(cp, "TPV.", 4) != 0) continue;
115 cp += 4;
116
117 sscanf(cp, "%zu", &index);
118 pv1[index]=disseq->dp[i].value.f;
119
120 /* For a check
121 printf("PV1_%ld\t%.10f\n", index, pv1[index]);
122 */
123 }
124
125 /* For axis2. */
126 else if (keyp->j == 2)
127 {
128 if ( keyp->field[0] == '\0' ) continue;
129 cp = strchr(keyp->field, '.') + 1;
130 if (strncmp(cp, "TPV.", 4) != 0) continue;
131 cp += 4;
132
133 sscanf(cp, "%zu", &index);
134
135 pv2[index]=disseq->dp[i].value.f;
136
137 /* For a check
138 printf("PV2_%ld\t%.10f\n", index, pv2[index]);
139 */
140 }
141 else
142 {
143 error(EXIT_FAILURE, 0, "%s: No such axis present!", __func__);
144 break;
145 }
146 }
147
148 /* To check a full array:
149 for(i = 0; i < 40; ++i)
150 printf("PV2_%ld\t%.10f\n", i, pv2[i]);
151 */
152
153 }
154
155
156
157
158
159
160 /* Extract the required sip parameters from wcs. */
161 static void
wcsdistortion_get_sipparams(struct wcsprm * wcs,double cd[2][2],double a_coeff[5][5],double b_coeff[5][5])162 wcsdistortion_get_sipparams(struct wcsprm *wcs, double cd[2][2],
163 double a_coeff[5][5], double b_coeff[5][5])
164 {
165 const char *cp;
166 double *temp_cd=NULL;
167 struct dpkey *keyp=NULL;
168 struct disprm *dispre=NULL;
169 size_t i=0, j=0, m=0, n=0, k=0;
170
171 /* Make sure a WCS structure is actually given. Note that 'wcs->lin' is
172 not a pointer, but the full structure. So it is always present. */
173 if(wcs==NULL)
174 error(EXIT_FAILURE, 0, "%s: input WCS structure is NULL", __func__);
175 if(wcs->lin.dispre == NULL)
176 error(EXIT_FAILURE, 0, "%s: input WCS structure's 'lin.dispre' is NULL",
177 __func__);
178
179 /* For easy reading. */
180 dispre=wcs->lin.dispre;
181 keyp=dispre->dp;
182
183 /* Fill the 2-times allocated CD array (cd[][]). Note that the required
184 CD matix is extracted using the `gal_wcs_wrap_matrix` as a single
185 allocated array (temp_cd[]), that is then used to fill cd[][]. */
186 temp_cd=gal_wcs_warp_matrix(wcs);
187 for(i=0,k=0; i<2; ++i)
188 for(j=0; j<2; ++j)
189 {
190 /* If a value is present store it, else store 0.*/
191 if(temp_cd[k] != 0) cd[i][j]=temp_cd[k++];
192 else cd[i][j]=0;
193
194 /* For a check:
195 printf("CD%ld_%ld\t%.10E\n", i+1, j+1, cd[i][j]);
196 */
197 }
198
199 /* Extract the SIP values from the strings and numbers inside the
200 wcsprm. */
201 for(i=0; i<dispre->ndp; ++i, ++keyp)
202 {
203 /* For axis1. */
204 if (keyp->j == 1)
205 {
206 if ( keyp->field[0] == '\0' ) continue;
207 cp = strchr(keyp->field, '.') + 1;
208 if (strncmp(cp, "SIP.FWD.", 8) != 0) continue;
209 cp += 8;
210
211 sscanf(cp, "%zu_%zu", &m, &n);
212 a_coeff[m][n]=dispre->dp[i].value.f;
213
214 /*For a check.
215 printf("A_%ld_%ld =\t %.10E\n", m, n, dispre->dp[i].value.f);
216 */
217 }
218
219 /* For axis2. */
220 else if (keyp->j == 2)
221 {
222 if ( keyp->field[0] == '\0' ) continue;
223 cp = strchr(keyp->field, '.') + 1;
224 if (strncmp(cp, "SIP.FWD.", 8) != 0) continue;
225 cp += 8;
226
227 sscanf(cp, "%zu_%zu", &m, &n);
228
229 b_coeff[m][n]=dispre->dp[i].value.f;
230
231 /*For a check.
232 printf("B_%ld_%ld =\t %.10E\n", m, n, dispre->dp[i].value.f);
233 */
234 }
235 else
236 error(EXIT_FAILURE, 0, "%s: No such axis present!", __func__);
237 }
238 }
239
240
241
242
243
244
245 /* Extract sip coefficients from the polynomial equation. */
246 static void
wcsdistortion_get_sipcoeff(struct wcsprm * wcs,size_t * a_order,size_t * b_order,double a_coeff[5][5],double b_coeff[5][5])247 wcsdistortion_get_sipcoeff(struct wcsprm *wcs, size_t *a_order,
248 size_t *b_order, double a_coeff[5][5],
249 double b_coeff[5][5])
250 {
251 size_t m, n;
252 double val=0;
253 size_t i, j, k;
254 double cd[2][2];
255 double tpvu[8][8], tpvv[8][8];
256
257 /* Initialise the 2d matrices. */
258 for(i=0; i<2; ++i) for(j=0; j<2; ++j) cd[i][j]=0;
259 for(i=0; i<8; ++i) for(j=0; j<8; ++j) { tpvu[i][j]=0; tpvv[i][j]=0; }
260
261 /* Calculate the TPV elements. */
262 wcsdistortion_calc_tpveq(wcs, cd, tpvu, tpvv);
263
264 for(i=0,k=0; i<2; ++i)
265 for(j=0; j<2; ++j)
266 {
267 /* If a value is present store it, else store 0.*/
268 if((wcs->cd)[i] != 0) cd[i][j]=(wcs->cd)[k++];
269 else cd[i][j]=0;
270
271 /* For a check:
272 printf("CD%ld_%ld\t%.10lf\n", i+1, j+1, cd[i][j]);
273 */
274 }
275
276 for(m=0; m<=4; ++m)
277 for(n=0; n<=4; ++n)
278 {
279 /*For axis = 1*/
280 val=wcsdistortion_calcsip(1, m, n, tpvu, tpvv);
281 if(val != 0)
282 {
283 a_coeff[m][n]=val;
284 *a_order=wcsdistortion_max(*a_order, wcsdistortion_max(m,n));
285 }
286 else
287 a_coeff[m][n]=0;
288
289 /*For axis = 2*/
290 val=wcsdistortion_calcsip(2, m, n, tpvu, tpvv);
291 if(val != 0)
292 {
293 b_coeff[m][n]=val;
294 *b_order=wcsdistortion_max(*b_order, wcsdistortion_max(m,n));
295 }
296 else
297 b_coeff[m][n]=0;
298 }
299
300 /*For a check.
301 for(m=0;m<=4;++m)
302 for(n=0;n<=4;++n)
303 printf("A_%ld_%ld = %.10E \t B_%ld_%ld = %.10E\n",
304 m, n, a_coeff[m][n], m, n, b_coeff[m][n]);
305 */
306 }
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327 /**************************************************************/
328 /********** Intermidiate Equations ************/
329 /**************************************************************/
330 /* Intermidiate equations formed during pv->sip conversions. */
331 static void
wcsdistortion_intermidate_tpveq(double cd[2][2],double * pv1,double * pv2,double k[5][5],double l[5][5])332 wcsdistortion_intermidate_tpveq(double cd[2][2], double *pv1, double *pv2,
333 double k[5][5], double l[5][5])
334 {
335 /*Intermidate polynomials, k[i][j] and l[i][j].
336 See Appendix A of the SPIE proceedings paper at
337 http://web.ipac.caltech.edu/staff/shupe/reprints/SIP_to_PV_SPIE2012.pdf
338
339 The work described in that paper is extended to 7th order.
340
341 References and citations:
342 "More flexibility in representing geometric distortion in
343 astronomical images,"
344 Shupe, David L.; Laher, Russ R.; Storrie-Lombardi, Lisa; Surace, Jason;
345 Grillmair, Carl; Levitan, David; Sesar, Branimir, 2012, in Software and
346 Cyberinfrastructure for Astronomy II.
347 Proceedings of the SPIE, Volume 8451, article id. 84511M.
348 */
349
350 k[0][0] = pv1[0];
351 l[0][0] = pv2[0];
352
353 k[0][1] = ( cd[0][1] * pv1[1]
354 + cd[1][1] * pv1[2] );
355 l[0][1] = ( cd[0][1] * pv2[2]
356 + cd[1][1] * pv2[1] );
357
358
359 k[1][0]= ( cd[0][0] * pv1[1]
360 + cd[1][0] * pv1[2] );
361 l[1][0]= ( cd[0][0] * pv2[2]
362 + cd[1][0] * pv2[1] );
363
364
365 k[0][2] = ( cd[0][1] * cd[0][1] * pv1[4]
366 + cd[0][1] * cd[1][1] * pv1[5]
367 + cd[1][1] * cd[1][1] * pv1[6] );
368 l[0][2] = ( cd[0][1] * cd[0][1] * pv2[6]
369 + cd[0][1] * cd[1][1] * pv2[5]
370 + cd[1][1] * cd[1][1] * pv2[4] );
371
372
373 k[1][1] = ( 2*cd[0][0] * cd[0][1] * pv1[4]
374 + cd[0][0] * cd[1][1] * pv1[5]
375 + cd[0][1] * cd[1][0] * pv1[5]
376 + 2*cd[1][0] * cd[1][1] * pv1[6] );
377 l[1][1] = ( 2*cd[0][0] * cd[0][1] * pv2[6]
378 + cd[0][0] * cd[1][1] * pv2[5]
379 + cd[0][1] * cd[1][0] * pv2[5]
380 + 2*cd[1][0] * cd[1][1] * pv2[4] );
381
382
383 k[2][0] = ( cd[0][0] * cd[0][0] * pv1[4]
384 + cd[0][0] * cd[1][0] * pv1[5]
385 + cd[1][0] * cd[1][0] * pv1[6] );
386 l[2][0] = ( cd[0][0] * cd[0][0] * pv2[6]
387 + cd[0][0] * cd[1][0] * pv2[5]
388 + cd[1][0] * cd[1][0] * pv2[4] );
389
390
391 k[0][3] = ( cd[0][1] * cd[0][1] * cd[0][1] * pv1[7]
392 + cd[0][1] * cd[0][1] * cd[1][1] * pv1[8]
393 + cd[0][1] * cd[1][1] * cd[1][1] * pv1[9]
394 + cd[1][1] * cd[1][1] * cd[1][1] * pv1[10] );
395 l[0][3] = ( cd[0][1] * cd[0][1] * cd[0][1] * pv2[10]
396 + cd[0][1] * cd[0][1] * cd[1][1] * pv2[9]
397 + cd[0][1] * cd[1][1] * cd[1][1] * pv2[8]
398 + cd[1][1] * cd[1][1] * cd[1][1] * pv2[7] );
399
400
401 k[1][2] = ( 3*cd[0][0] * cd[0][1] * cd[0][1] * pv1[7]
402 + 2*cd[0][0] * cd[0][1] * cd[1][1] * pv1[8]
403 + cd[0][0] * cd[1][1] * cd[1][1] * pv1[9]
404 + cd[0][1] * cd[0][1] * cd[1][0] * pv1[8]
405 + 2*cd[0][1] * cd[1][0] * cd[1][1] * pv1[9]
406 + 3*cd[1][0] * cd[1][1] * cd[1][1] * pv1[10] );
407 l[1][2] = ( 3*cd[0][0] * cd[0][1] * cd[0][1] * pv2[10]
408 + 2*cd[0][0] * cd[0][1] * cd[1][1] * pv2[9]
409 + cd[0][0] * cd[1][1] * cd[1][1] * pv2[8]
410 + cd[0][1] * cd[0][1] * cd[1][0] * pv2[9]
411 + 2*cd[0][1] * cd[1][0] * cd[1][1] * pv2[8]
412 + 3*cd[1][0] * cd[1][1] * cd[1][1] * pv2[7] );
413
414
415 k[2][1] = ( 3*cd[0][0] * cd[0][0] * cd[0][1] * pv1[7]
416 + 2*cd[0][0] * cd[0][1] * cd[1][0] * pv1[8]
417 + cd[0][0] * cd[0][0] * cd[1][1] * pv1[8]
418 + cd[0][1] * cd[1][0] * cd[1][0] * pv1[9]
419 + 2*cd[0][0] * cd[1][0] * cd[1][1] * pv1[9]
420 + 3*cd[1][0] * cd[1][0] * cd[1][1] * pv1[10] );
421 l[2][1] = ( 3*cd[0][0] * cd[0][0] * cd[0][1] * pv2[10]
422 + 2*cd[0][0] * cd[0][1] * cd[1][0] * pv2[9]
423 + cd[0][0] * cd[0][0] * cd[1][1] * pv2[9]
424 + cd[0][1] * cd[1][0] * cd[1][0] * pv2[8]
425 + 2*cd[0][0] * cd[1][0] * cd[1][1] * pv2[8]
426 + 3*cd[1][0] * cd[1][0] * cd[1][1] * pv2[7] );
427
428
429 k[3][0] = ( cd[0][0] * cd[0][0] * cd[0][0] * pv1[7]
430 + cd[0][0] * cd[0][0] * cd[1][0] * pv1[8]
431 + cd[0][0] * cd[1][0] * cd[1][0] * pv1[9]
432 + cd[1][0] * cd[1][0] * cd[1][0] * pv1[10] );
433 l[3][0] = ( cd[0][0] * cd[0][0] * cd[0][0] * pv2[10]
434 + cd[0][0] * cd[0][0] * cd[1][0] * pv2[9]
435 + cd[0][0] * cd[1][0] * cd[1][0] * pv2[8]
436 + cd[1][0] * cd[1][0] * cd[1][0] * pv2[7] );
437
438
439 k[0][4] = ( cd[0][1] * cd[0][1] * cd[0][1] * cd[0][1] * pv1[12]
440 + cd[0][1] * cd[0][1] * cd[1][1] * cd[1][1] * pv1[13]
441 + cd[0][1] * cd[0][1] * cd[1][1] * cd[1][1] * pv1[14]
442 + cd[0][1] * cd[1][1] * cd[1][1] * cd[1][1] * pv1[15]
443 + cd[1][1] * cd[1][1] * cd[1][1] * cd[1][1] * pv1[16] );
444 l[0][4] = ( cd[0][1] * cd[0][1] * cd[0][1] * cd[0][1] * pv2[16]
445 + cd[0][1] * cd[0][1] * cd[1][1] * cd[1][1] * pv2[15]
446 + cd[0][1] * cd[0][1] * cd[1][1] * cd[1][1] * pv2[14]
447 + cd[0][1] * cd[1][1] * cd[1][1] * cd[1][1] * pv2[13]
448 + cd[1][1] * cd[1][1] * cd[1][1] * cd[1][1] * pv2[12] );
449
450
451 k[1][3] = ( 4*cd[0][0] * cd[0][1] * cd[0][1] * cd[0][1] * pv1[12]
452 + 3*cd[0][0] * cd[0][1] * cd[0][1] * cd[1][1] * pv1[13]
453 + 2*cd[0][0] * cd[0][1] * cd[1][1] * cd[1][1] * pv1[14]
454 + cd[0][0] * cd[1][1] * cd[1][1] * cd[1][1] * pv1[15]
455 + cd[0][1] * cd[0][1] * cd[0][1] * cd[1][0] * pv1[13]
456 + 2*cd[0][1] * cd[0][1] * cd[1][0] * cd[1][1] * pv1[14]
457 + 3*cd[0][1] * cd[1][0] * cd[1][1] * cd[1][1] * pv1[15]
458 + 4*cd[1][0] * cd[1][1] * cd[1][1] * cd[1][1] * pv1[16] );
459 l[1][3] = ( 4*cd[0][0] * cd[0][1] * cd[0][1] * cd[0][1] * pv2[16]
460 + 3*cd[0][0] * cd[0][1] * cd[0][1] * cd[1][1] * pv2[15]
461 + 2*cd[0][0] * cd[0][1] * cd[1][1] * cd[1][1] * pv2[14]
462 + cd[0][0] * cd[1][1] * cd[1][1] * cd[1][1] * pv2[13]
463 + cd[0][1] * cd[0][1] * cd[0][1] * cd[1][0] * pv2[15]
464 + 2*cd[0][1] * cd[0][1] * cd[1][0] * cd[1][1] * pv2[14]
465 + 3*cd[0][1] * cd[1][0] * cd[1][1] * cd[1][1] * pv2[13]
466 + 4*cd[1][0] * cd[1][1] * cd[1][1] * cd[1][1] * pv2[12] );
467
468
469 k[2][2] = ( 6*cd[0][0] * cd[0][0] * cd[0][1] * cd[0][1] * pv1[12]
470 + 3*cd[0][0] * cd[0][0] * cd[0][1] * cd[1][1] * pv1[13]
471 + cd[0][0] * cd[0][0] * cd[1][1] * cd[1][1] * pv1[14]
472 + 3*cd[0][0] * cd[0][1] * cd[0][1] * cd[1][0] * pv1[13]
473 + 4*cd[0][0] * cd[0][1] * cd[1][0] * cd[1][1] * pv1[14]
474 + 3*cd[0][0] * cd[1][0] * cd[1][1] * cd[1][1] * pv1[15]
475 + cd[0][1] * cd[0][1] * cd[1][0] * cd[1][0] * pv1[14]
476 + 3*cd[0][1] * cd[1][0] * cd[1][0] * cd[1][1] * pv1[15]
477 + 6*cd[1][0] * cd[1][0] * cd[1][1] * cd[1][1] * pv1[16] );
478 l[2][2] = ( 6*cd[0][0] * cd[0][0] * cd[0][1] * cd[0][1] * pv2[16]
479 + 3*cd[0][0] * cd[0][0] * cd[0][1] * cd[1][1] * pv2[15]
480 + cd[0][0] * cd[0][0] * cd[1][1] * cd[1][1] * pv2[14]
481 + 3*cd[0][0] * cd[0][1] * cd[0][1] * cd[1][0] * pv2[15]
482 + 4*cd[0][0] * cd[0][1] * cd[1][0] * cd[1][1] * pv2[14]
483 + 3*cd[0][0] * cd[1][0] * cd[1][1] * cd[1][1] * pv2[15]
484 + cd[0][1] * cd[0][1] * cd[1][0] * cd[1][0] * pv2[14]
485 + 3*cd[0][1] * cd[1][0] * cd[1][0] * cd[1][1] * pv2[13]
486 + 6*cd[1][0] * cd[1][0] * cd[1][1] * cd[1][1] * pv2[12] );
487
488
489 k[3][1] = ( 4*cd[0][0] * cd[0][0] * cd[0][0] * cd[0][1] * pv1[12]
490 + cd[0][0] * cd[0][0] * cd[0][0] * cd[1][1] * pv1[13]
491 + 3*cd[0][0] * cd[0][0] * cd[0][1] * cd[1][0] * pv1[13]
492 + 2*cd[0][0] * cd[0][0] * cd[1][0] * cd[1][1] * pv1[14]
493 + 2*cd[0][0] * cd[0][1] * cd[1][0] * cd[1][0] * pv1[14]
494 + 3*cd[0][0] * cd[1][0] * cd[1][0] * cd[1][1] * pv1[15]
495 + cd[0][1] * cd[1][0] * cd[1][0] * cd[1][0] * pv1[15]
496 + 4*cd[1][0] * cd[1][0] * cd[1][0] * cd[1][1] * pv1[16] );
497 l[3][1] = ( 4*cd[0][0] * cd[0][0] * cd[0][0] * cd[0][1] * pv2[16]
498 + cd[0][0] * cd[0][0] * cd[0][0] * cd[1][1] * pv2[15]
499 + 3*cd[0][0] * cd[0][0] * cd[0][1] * cd[1][0] * pv2[15]
500 + 2*cd[0][0] * cd[0][0] * cd[1][0] * cd[1][1] * pv2[14]
501 + 2*cd[0][0] * cd[0][1] * cd[1][0] * cd[1][0] * pv2[14]
502 + 3*cd[0][0] * cd[1][0] * cd[1][0] * cd[1][1] * pv2[13]
503 + cd[0][1] * cd[1][0] * cd[1][0] * cd[1][0] * pv2[13]
504 + 4*cd[1][0] * cd[1][0] * cd[1][0] * cd[1][1] * pv2[12] );
505
506
507 k[4][0] = ( cd[0][0] * cd[0][0] * cd[0][0] * cd[0][0] * pv1[12]
508 + cd[0][0] * cd[0][0] * cd[0][0] * cd[1][0] * pv1[13]
509 + cd[0][0] * cd[0][0] * cd[1][0] * cd[1][0] * pv1[14]
510 + cd[0][0] * cd[1][0] * cd[1][0] * cd[1][0] * pv1[15]
511 + cd[1][0] * cd[1][0] * cd[1][0] * cd[1][0] * pv1[16] );
512 l[4][0] = ( cd[0][0] * cd[0][0] * cd[0][0] * cd[0][0] * pv2[16]
513 + cd[0][0] * cd[0][0] * cd[0][0] * cd[1][0] * pv2[15]
514 + cd[0][0] * cd[0][0] * cd[1][0] * cd[1][0] * pv2[14]
515 + cd[0][0] * cd[1][0] * cd[1][0] * cd[1][0] * pv2[13]
516 + cd[1][0] * cd[1][0] * cd[1][0] * cd[1][0] * pv2[12] );
517
518
519 /* For a check:
520 {
521 size_t i, j;
522 for(i=0; i<=4; ++i)
523 for(j=0;j<=4;++j)
524 {
525 printf("k%ld_%ld \t %.10E\n", i, j, k[i][j]);
526 printf("l%ld_%ld \t %.10E\n\n", i, j, l[i][j]);
527 }
528 }
529 */
530 }
531
532
533
534
535
536
537 /* Intermidiate equations formed during sip->pv conversions. */
538 static void
wcsdistortion_intermidate_sipeq(double cd[2][2],double cd_inv[2][2],double a_coeff[5][5],double b_coeff[5][5],double * pv1,double * pv2)539 wcsdistortion_intermidate_sipeq(double cd[2][2], double cd_inv[2][2],
540 double a_coeff[5][5], double b_coeff[5][5],
541 double *pv1, double *pv2)
542 {
543 /*Intemidiate equations which give the value of PVi_j
544 excluding radial terms PV[i,3], PV[i,11]
545
546 For the intermidiate PVi_j calculations, the Eqs. 1 and 2 from
547 section 2 of the SPIE proceedings paper at
548 http://web.ipac.caltech.edu/staff/shupe/reprints/SIP_to_PV_SPIE2012.pdf
549 are used which were further modified as shown in
550 https://github.com/stargaser/sip_tpv/blob/master/sip_tpv/pvsiputils.py
551 to generate these equtions.
552
553 The exact script used for generation of these equations are given here:
554 https://gitlab.com/sachinkumarsingh092/gnuastro-test-files/-/blob/
555 master/scripts/equations.py
556
557 References and citations:
558 "More flexibility in representing geometric distortion in
559 astronomical images,"
560 Shupe, David L.; Laher, Russ R.; Storrie-Lombardi, Lisa; Surace, Jason;
561 Grillmair, Carl; Levitan, David; Sesar, Branimir, 2012, in Software and
562 Cyberinfrastructure for Astronomy II.
563 Proceedings of the SPIE, Volume 8451, article id. 84511M.
564 */
565
566 /* pvi_0 */
567 pv1[0]=a_coeff[0][0]*cd[0][0] + b_coeff[0][0]*cd[0][1];
568 pv2[0]=a_coeff[0][0]*cd[1][0] + b_coeff[0][0]*cd[1][1];
569
570
571 /* pvi_1 */
572 pv1[1] = ( a_coeff[0][1] * cd[0][0] * cd_inv[1][0]
573 + a_coeff[1][0] * cd[0][0] * cd_inv[0][0]
574 + b_coeff[0][1] * cd[0][1] * cd_inv[1][0]
575 + b_coeff[1][0] * cd[0][1] * cd_inv[0][0]
576 + cd[0][0] * cd_inv[0][0]
577 + cd[0][1] * cd_inv[1][0] );
578
579 pv2[1] = ( a_coeff[0][1] * cd[1][0] * cd_inv[1][1]
580 + a_coeff[1][0] * cd[1][0] * cd_inv[0][1]
581 + b_coeff[0][1] * cd[1][1] * cd_inv[1][1]
582 + b_coeff[1][0] * cd[1][1] * cd_inv[0][1]
583 + cd[1][0] * cd_inv[0][1]
584 + cd[1][1] * cd_inv[1][1] );
585
586
587 /* pvi_2 */
588 pv1[2] = ( a_coeff[0][1] * cd[0][0] * cd_inv[1][1]
589 + a_coeff[1][0] * cd[0][0] * cd_inv[0][1]
590 + b_coeff[0][1] * cd[0][1] * cd_inv[1][1]
591 + b_coeff[1][0] * cd[0][1] * cd_inv[0][1]
592 + cd[0][0] * cd_inv[0][1]
593 + cd[0][1] * cd_inv[1][1] );
594
595 pv2[2] = ( a_coeff[0][1] * cd[1][0] * cd_inv[1][0]
596 + a_coeff[1][0] * cd[1][0] * cd_inv[0][0]
597 + b_coeff[0][1] * cd[1][1] * cd_inv[1][0]
598 + b_coeff[1][0] * cd[1][1] * cd_inv[0][0]
599 + cd[1][0] * cd_inv[0][0]
600 + cd[1][1] * cd_inv[1][0] );
601
602
603 /* pvi_4 */
604 pv1[4] = ( a_coeff[0][2] * cd[0][0] * cd_inv[1][0] * cd_inv[1][0]
605 + a_coeff[1][1] * cd[0][0] * cd_inv[0][0] * cd_inv[1][0]
606 + a_coeff[2][0] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0]
607 + b_coeff[0][2] * cd[0][1] * cd_inv[1][0] * cd_inv[1][0]
608 + b_coeff[1][1] * cd[0][1] * cd_inv[0][0] * cd_inv[1][0]
609 + b_coeff[2][0] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] );
610
611 pv2[4] = ( a_coeff[0][2] * cd[1][0] * cd_inv[1][1] * cd_inv[1][1]
612 + a_coeff[1][1] * cd[1][0] * cd_inv[0][1] * cd_inv[1][1]
613 + a_coeff[2][0] * cd[1][0] * cd_inv[0][1] * cd_inv[0][1]
614 + b_coeff[0][2] * cd[1][1] * cd_inv[1][1] * cd_inv[1][1]
615 + b_coeff[1][1] * cd[1][1] * cd_inv[0][1] * cd_inv[1][1]
616 + b_coeff[2][0] * cd[1][1] * cd_inv[0][1] * cd_inv[0][1] );
617
618
619 /* pvi_5 */
620 pv1[5] = ( 2*a_coeff[0][2] * cd[0][0] * cd_inv[1][0] * cd_inv[1][1]
621 + a_coeff[1][1] * cd[0][0] * cd_inv[0][0] * cd_inv[1][1]
622 + a_coeff[1][1] * cd[0][0] * cd_inv[0][1] * cd_inv[1][0]
623 + 2*a_coeff[2][0] * cd[0][0] * cd_inv[0][0] * cd_inv[0][1]
624 + 2*b_coeff[0][2] * cd[0][1] * cd_inv[1][0] * cd_inv[1][1]
625 + b_coeff[1][1] * cd[0][1] * cd_inv[0][0] * cd_inv[1][1]
626 + b_coeff[1][1] * cd[0][1] * cd_inv[0][1] * cd_inv[1][0]
627 + 2*b_coeff[2][0] * cd[0][1] * cd_inv[0][0] * cd_inv[0][1] );
628
629 pv2[5] = ( 2*a_coeff[0][2] * cd[1][0] * cd_inv[1][0] * cd_inv[1][1]
630 + a_coeff[1][1] * cd[1][0] * cd_inv[0][0] * cd_inv[1][1]
631 + a_coeff[1][1] * cd[1][0] * cd_inv[0][1] * cd_inv[1][0]
632 + 2*a_coeff[2][0] * cd[1][0] * cd_inv[0][0] * cd_inv[0][1]
633 + 2*b_coeff[0][2] * cd[1][1] * cd_inv[1][0] * cd_inv[1][1]
634 + b_coeff[1][1] * cd[1][1] * cd_inv[0][0] * cd_inv[1][1]
635 + b_coeff[1][1] * cd[1][1] * cd_inv[0][1] * cd_inv[1][0]
636 + 2*b_coeff[2][0] * cd[1][1] * cd_inv[0][0] * cd_inv[0][1] );
637
638
639 /* pvi_6 */
640 pv1[6]= ( a_coeff[0][2] * cd[0][0] * cd_inv[1][1] * cd_inv[1][1]
641 + a_coeff[1][1] * cd[0][0] * cd_inv[0][1] * cd_inv[1][1]
642 + a_coeff[2][0] * cd[0][0] * cd_inv[0][1] * cd_inv[0][1]
643 + b_coeff[0][2] * cd[0][1] * cd_inv[1][1] * cd_inv[1][1]
644 + b_coeff[1][1] * cd[0][1] * cd_inv[0][1] * cd_inv[1][1]
645 + b_coeff[2][0] * cd[0][1] * cd_inv[0][1] * cd_inv[0][1] );
646
647 pv2[6]= ( a_coeff[0][2] * cd[1][0] * cd_inv[1][0] * cd_inv[1][0]
648 + a_coeff[1][1] * cd[1][0] * cd_inv[0][0] * cd_inv[1][0]
649 + a_coeff[2][0] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0]
650 + b_coeff[0][2] * cd[1][1] * cd_inv[1][0] * cd_inv[1][0]
651 + b_coeff[1][1] * cd[1][1] * cd_inv[0][0] * cd_inv[1][0]
652 + b_coeff[2][0] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] );
653
654
655 /* pvi_7 */
656 pv1[7] = ( a_coeff[0][3] * cd[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
657 + a_coeff[1][2] * cd[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0]
658 + a_coeff[2][1] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0]
659 + a_coeff[3][0] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0]
660 + b_coeff[0][3] * cd[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
661 + b_coeff[1][2] * cd[0][1] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0]
662 + b_coeff[2][1] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0]
663 + b_coeff[3][0] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] );
664
665 pv2[7] = ( a_coeff[0][3] * cd[1][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
666 + a_coeff[1][2] * cd[1][0] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
667 + a_coeff[2][1] * cd[1][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
668 + a_coeff[3][0] * cd[1][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1]
669 + b_coeff[0][3] * cd[1][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
670 + b_coeff[1][2] * cd[1][1] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
671 + b_coeff[2][1] * cd[1][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
672 + b_coeff[3][0] * cd[1][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] );
673
674
675 /* pvi_8 */
676 pv1[8] = ( 3*a_coeff[0][3] * cd[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
677 + 2*a_coeff[1][2] * cd[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1]
678 + a_coeff[1][2] * cd[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
679 + a_coeff[2][1] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1]
680 + 2*a_coeff[2][1] * cd[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0]
681 + 3*a_coeff[3][0] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1]
682 + 3*b_coeff[0][3] * cd[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
683 + 2*b_coeff[1][2] * cd[0][1] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1]
684 + b_coeff[1][2] * cd[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
685 + b_coeff[2][1] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1]
686 + 2*b_coeff[2][1] * cd[0][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0]
687 + 3*b_coeff[3][0] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] );
688
689 pv2[8] = ( 3*a_coeff[0][3] * cd[1][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
690 + a_coeff[1][2] * cd[1][0] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1]
691 + 2*a_coeff[1][2] * cd[1][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
692 + 2*a_coeff[2][1] * cd[1][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1]
693 + a_coeff[2][1] * cd[1][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
694 + 3*a_coeff[3][0] * cd[1][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1]
695 + 3*b_coeff[0][3] * cd[1][1] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
696 + b_coeff[1][2] * cd[1][1] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1]
697 + 2*b_coeff[1][2] * cd[1][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
698 + 2*b_coeff[2][1] * cd[1][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1]
699 + b_coeff[2][1] * cd[1][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
700 + 3*b_coeff[3][0] * cd[1][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] );
701
702
703 /* pvi_9 */
704 pv1[9] = ( 3*a_coeff[0][3] * cd[0][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
705 + a_coeff[1][2] * cd[0][0] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1]
706 + 2*a_coeff[1][2] * cd[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
707 + 2*a_coeff[2][1] * cd[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1]
708 + a_coeff[2][1] * cd[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
709 + 3*a_coeff[3][0] * cd[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1]
710 + 3*b_coeff[0][3] * cd[0][1] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
711 + b_coeff[1][2] * cd[0][1] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1]
712 + 2*b_coeff[1][2] * cd[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
713 + 2*b_coeff[2][1] * cd[0][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1]
714 + b_coeff[2][1] * cd[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
715 + 3*b_coeff[3][0] * cd[0][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] );
716
717 pv2[9] = ( 3*a_coeff[0][3] * cd[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
718 + 2*a_coeff[1][2] * cd[1][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1]
719 + a_coeff[1][2] * cd[1][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
720 + a_coeff[2][1] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1]
721 + 2*a_coeff[2][1] * cd[1][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0]
722 + 3*a_coeff[3][0] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1]
723 + 3*b_coeff[0][3] * cd[1][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
724 + 2*b_coeff[1][2] * cd[1][1] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1]
725 + b_coeff[1][2] * cd[1][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
726 + b_coeff[2][1] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1]
727 + 2*b_coeff[2][1] * cd[1][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0]
728 + 3*b_coeff[3][0] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] );
729
730
731 /* pvi_10 */
732 pv1[10] = ( a_coeff[0][3] * cd[0][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
733 + a_coeff[1][2] * cd[0][0] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
734 + a_coeff[2][1] * cd[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
735 + a_coeff[3][0] * cd[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1]
736 + b_coeff[0][3] * cd[0][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
737 + b_coeff[1][2] * cd[0][1] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
738 + b_coeff[2][1] * cd[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
739 + b_coeff[3][0] * cd[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] );
740
741
742 pv2[10] = ( a_coeff[0][3] * cd[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
743 + a_coeff[1][2] * cd[1][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0]
744 + a_coeff[2][1] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0]
745 + a_coeff[3][0] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0]
746 + b_coeff[0][3] * cd[1][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
747 + b_coeff[1][2] * cd[1][1] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0]
748 + b_coeff[2][1] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0]
749 + b_coeff[3][0] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] );
750
751
752 /* pvi_12 */
753 pv1[12] = ( a_coeff[0][4] * cd[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
754 + a_coeff[1][3] * cd[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
755 + a_coeff[2][2] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0]
756 + a_coeff[3][1] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0]
757 + a_coeff[4][0] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0]
758 + b_coeff[0][4] * cd[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
759 + b_coeff[1][3] * cd[0][1] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
760 + b_coeff[2][2] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0]
761 + b_coeff[3][1] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0]
762 + b_coeff[4][0] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] );
763
764 pv2[12] = ( a_coeff[0][4] * cd[1][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
765 + a_coeff[1][3] * cd[1][0] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
766 + a_coeff[2][2] * cd[1][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
767 + a_coeff[3][1] * cd[1][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
768 + a_coeff[4][0] * cd[1][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1]
769 + b_coeff[0][4] * cd[1][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
770 + b_coeff[1][3] * cd[1][1] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
771 + b_coeff[2][2] * cd[1][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
772 + b_coeff[3][1] * cd[1][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
773 + b_coeff[4][0] * cd[1][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] );
774
775
776 /* pvi_13 */
777 pv1[13] = ( 4*a_coeff[0][4] * cd[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
778 + 3*a_coeff[1][3] * cd[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
779 + a_coeff[1][3] * cd[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
780 + 2*a_coeff[2][2] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1]
781 + 2*a_coeff[2][2] * cd[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
782 + a_coeff[3][1] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1]
783 + 3*a_coeff[3][1] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0]
784 + 4*a_coeff[4][0] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1]
785 + 4*b_coeff[0][4] * cd[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
786 + 3*b_coeff[1][3] * cd[0][1] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
787 + b_coeff[1][3] * cd[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
788 + 2*b_coeff[2][2] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1]
789 + 2*b_coeff[2][2] * cd[0][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
790 + b_coeff[3][1] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1]
791 + 3*b_coeff[3][1] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0]
792 + 4*b_coeff[4][0] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] );
793
794
795 pv2[13] = ( 4*a_coeff[0][4] * cd[1][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
796 + a_coeff[1][3] * cd[1][0] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
797 + 3*a_coeff[1][3] * cd[1][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
798 + 2*a_coeff[2][2] * cd[1][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
799 + 2*a_coeff[2][2] * cd[1][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
800 + 3*a_coeff[3][1] * cd[1][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
801 + a_coeff[3][1] * cd[1][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
802 + 4*a_coeff[4][0] * cd[1][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1]
803 + 4*b_coeff[0][4] * cd[1][1] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
804 + b_coeff[1][3] * cd[1][1] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
805 + 3*b_coeff[1][3] * cd[1][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
806 + 2*b_coeff[2][2] * cd[1][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
807 + 2*b_coeff[2][2] * cd[1][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
808 + 3*b_coeff[3][1] * cd[1][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
809 + b_coeff[3][1] * cd[1][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
810 + 4*b_coeff[4][0] * cd[1][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] );
811
812
813 /* pvi_14 */
814 pv1[14] = ( 6*a_coeff[0][4] * cd[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
815 + 3*a_coeff[1][3] * cd[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
816 + 3*a_coeff[1][3] * cd[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
817 + a_coeff[2][2] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1]
818 + 4*a_coeff[2][2] * cd[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
819 + a_coeff[2][2] * cd[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
820 + 3*a_coeff[3][1] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1]
821 + 3*a_coeff[3][1] * cd[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
822 + 6*a_coeff[4][0] * cd[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1]
823 + 6*b_coeff[0][4] * cd[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
824 + 3*b_coeff[1][3] * cd[0][1] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
825 + 3*b_coeff[1][3] * cd[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
826 + b_coeff[2][2] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1]
827 + 4*b_coeff[2][2] * cd[0][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
828 + b_coeff[2][2] * cd[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
829 + 3*b_coeff[3][1] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1]
830 + 3*b_coeff[3][1] * cd[0][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
831 + 6*b_coeff[4][0] * cd[0][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] );
832
833 pv2[14] = ( 6*a_coeff[0][4] * cd[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
834 + 3*a_coeff[1][3] * cd[1][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
835 + 3*a_coeff[1][3] * cd[1][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
836 + a_coeff[2][2] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1]
837 + 4*a_coeff[2][2] * cd[1][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
838 + a_coeff[2][2] * cd[1][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
839 + 3*a_coeff[3][1] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1]
840 + 3*a_coeff[3][1] * cd[1][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
841 + 6*a_coeff[4][0] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1]
842 + 6*b_coeff[0][4] * cd[1][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
843 + 3*b_coeff[1][3] * cd[1][1] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
844 + 3*b_coeff[1][3] * cd[1][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
845 + b_coeff[2][2] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1]
846 + 4*b_coeff[2][2] * cd[1][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
847 + b_coeff[2][2] * cd[1][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
848 + 3*b_coeff[3][1] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1]
849 + 3*b_coeff[3][1] * cd[1][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
850 + 6*b_coeff[4][0] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] );
851
852
853 /* pvi_15 */
854 pv1[15] = ( 4*a_coeff[0][4] * cd[0][0] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
855 + a_coeff[1][3] * cd[0][0] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
856 + 3*a_coeff[1][3] * cd[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
857 + 2*a_coeff[2][2] * cd[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
858 + 2*a_coeff[2][2] * cd[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
859 + 3*a_coeff[3][1] * cd[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
860 + a_coeff[3][1] * cd[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
861 + 4*a_coeff[4][0] * cd[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1]
862 + 4*b_coeff[0][4] * cd[0][1] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
863 + b_coeff[1][3] * cd[0][1] * cd_inv[0][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
864 + 3*b_coeff[1][3] * cd[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1] * cd_inv[1][1]
865 + 2*b_coeff[2][2] * cd[0][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
866 + 2*b_coeff[2][2] * cd[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][1]
867 + 3*b_coeff[3][1] * cd[0][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
868 + b_coeff[3][1] * cd[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][0]
869 + 4*b_coeff[4][0] * cd[0][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] );
870
871 pv2[15] = ( 4*a_coeff[0][4] * cd[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
872 + 3*a_coeff[1][3] * cd[1][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
873 + a_coeff[1][3] * cd[1][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
874 + 2*a_coeff[2][2] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1]
875 + 2*a_coeff[2][2] * cd[1][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
876 + a_coeff[3][1] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1]
877 + 3*a_coeff[3][1] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0]
878 + 4*a_coeff[4][0] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1]
879 + 4*b_coeff[0][4] * cd[1][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
880 + 3*b_coeff[1][3] * cd[1][1] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][1]
881 + b_coeff[1][3] * cd[1][1] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
882 + 2*b_coeff[2][2] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][1]
883 + 2*b_coeff[2][2] * cd[1][1] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0] * cd_inv[1][0]
884 + b_coeff[3][1] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][1]
885 + 3*b_coeff[3][1] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] * cd_inv[1][0]
886 + 4*b_coeff[4][0] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][1] );
887
888
889 /* pvi_16 */
890 pv1[16] = ( a_coeff[0][4] * cd[0][0] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
891 + a_coeff[1][3] * cd[0][0] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
892 + a_coeff[2][2] * cd[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
893 + a_coeff[3][1] * cd[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
894 + a_coeff[4][0] * cd[0][0] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1]
895 + b_coeff[0][4] * cd[0][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
896 + b_coeff[1][3] * cd[0][1] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1] * cd_inv[1][1]
897 + b_coeff[2][2] * cd[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1] * cd_inv[1][1]
898 + b_coeff[3][1] * cd[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[1][1]
899 + b_coeff[4][0] * cd[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] * cd_inv[0][1] );
900
901 pv2[16] = ( a_coeff[0][4] * cd[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
902 + a_coeff[1][3] * cd[1][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
903 + a_coeff[2][2] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0]
904 + a_coeff[3][1] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0]
905 + a_coeff[4][0] * cd[1][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0]
906 + b_coeff[0][4] * cd[1][1] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
907 + b_coeff[1][3] * cd[1][1] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0] * cd_inv[1][0]
908 + b_coeff[2][2] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0] * cd_inv[1][0]
909 + b_coeff[3][1] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[1][0]
910 + b_coeff[4][0] * cd[1][1] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] * cd_inv[0][0] );
911
912
913 /* For a check:
914 {
915 size_t j;
916 for(j=0; j<=16 ;++j)
917 {
918 if (j == 3 || j == 11) continue;
919 printf("pv1_%ld \t %.12E\n", j, pv1[j]);
920 printf("pv2_%ld \t %.12E\n", j, pv2[j]);
921 }
922 }
923 */
924 }
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945 /**************************************************************/
946 /********** Calculations ************/
947 /**************************************************************/
948 /* Calcualte the final PV equation from intermidiate equations. */
949 static void
wcsdistortion_calc_tpveq(struct wcsprm * wcs,double cd[2][2],double tpvu[8][8],double tpvv[8][8])950 wcsdistortion_calc_tpveq(struct wcsprm *wcs, double cd[2][2],
951 double tpvu[8][8], double tpvv[8][8])
952 {
953 /* tpvu and tpvv are u-v distortion equations in TPV convention. */
954 size_t i=0,j=0;
955 double cd_inv[2][2];
956 double determinant=0;
957 double k[5][5], l[5][5];
958 double pv1[17]={0}, pv2[17]={0};
959
960 /* Initialise the 2d matrices. */
961 for(i=0; i<2; ++i) for(j=0; j<2; ++j) cd_inv[i][j]=0;
962 for(i=0; i<5; ++i) for(j=0; j<5; ++j) { k[i][j]=0; l[i][j]=0; }
963
964 /* Estimate the parameters. */
965 wcsdistortion_get_tpvparams(wcs, cd, pv1, pv2);
966 wcsdistortion_intermidate_tpveq(cd, pv1, pv2, k, l);
967
968 /* We will find matrix tpvu and tpvv by finding inverse of
969 CD matrix and multiplying with tpv* matrix.
970 For inverse of a 2x2 matrix we use the below trasformations:
971 inverse(|a b|) = 1 *|d -b|
972 |c d| |A| |-c a|
973 where |A| is the determinant of the matrix which is calculated by:
974 |A| = a*d-b*c.
975 */
976 determinant = cd[0][0]*cd[1][1] - cd[0][1]*cd[1][0];
977
978 /* Inverse matrix */
979 cd_inv[0][0]=cd[1][1]/determinant; /* a */
980 cd_inv[0][1]=(-1*cd[0][1])/determinant; /* b */
981 cd_inv[1][0]=(-1*cd[1][0])/determinant; /* c */
982 cd_inv[1][1]=cd[0][0]/determinant; /* d */
983
984 /*For a check.
985 printf("%.10lf\t%.10lf\t%.10lf\t%.10lf\n", cd_inv[0][0], cd_inv[0][1], \
986 cd_inv[1][0], cd_inv[1][1]);
987 printf("%.10lf\n", determinant);
988 */
989
990 /* For matrix tpvv and tpvu, we have to use the following
991 matrix equation:
992
993 |tpvu| = cd_inv*|k[i][j]|
994 |tpvv| |l[i][j]|
995 though intermidate polynomial equations have to be calculated prior
996 to multiplycation with cd_inv.
997 */
998
999 for(i=0; i<=4; ++i)
1000 for(j=0; j<=4; ++j)
1001 {
1002 tpvu[i][j]=cd_inv[0][0]*k[i][j]+cd_inv[0][1]*l[i][j];
1003 tpvv[i][j]=cd_inv[1][0]*k[i][j]+cd_inv[1][1]*l[i][j];
1004
1005 /*For a check:
1006 printf("%.10E, %.10E\n", tpvu[i][j], tpvv[i][j]);
1007 */
1008 }
1009 }
1010
1011
1012
1013
1014
1015 /* Calcualte the final SIP equation from intermidiate equations
1016 and extract pv coefficients from it. */
1017 static void
wcsdistortion_calc_sipeq(struct wcsprm * wcs,double cd[2][2],double * pv1,double * pv2)1018 wcsdistortion_calc_sipeq(struct wcsprm *wcs, double cd[2][2],
1019 double *pv1, double *pv2)
1020 {
1021 /* tpvu and tpvv are u-v distortion equations in TPV convention. */
1022 size_t i=0, j=0;
1023 double cd_inv[2][2];
1024 double determinant=0;
1025 double a_coeff[5][5], b_coeff[5][5];
1026
1027 /* Initialise the 2d matrices. */
1028 for(i=0;i<2;++i) for(j=0;j<2;++j) cd_inv[i][j]=0;
1029 for(i=0;i<5;++i) for(j=0;j<5;++j) {a_coeff[i][j]=0; b_coeff[i][j]=0;}
1030
1031 /* We will find matrix tpvu and tpvv by finding inverse of
1032 CD matrix and multiplying with tpv* matrix.
1033 For inverse of a 2x2 matrix we use the below trasformations:
1034 inverse(|a b|) = 1 *|d -b|
1035 |c d| |A| |-c a|
1036 where |A| is the determinant of the matrix which is calculated by:
1037 |A| = a*d-b*c.
1038 */
1039 wcsdistortion_get_sipparams(wcs, cd, a_coeff, b_coeff);
1040 determinant = cd[0][0]*cd[1][1] - cd[0][1]*cd[1][0];
1041
1042 /* Inverse matrix */
1043 cd_inv[0][0]=cd[1][1]/determinant; /* a */
1044 cd_inv[0][1]=(-1*cd[0][1])/determinant; /* b */
1045 cd_inv[1][0]=(-1*cd[1][0])/determinant; /* c */
1046 cd_inv[1][1]=cd[0][0]/determinant; /* d */
1047
1048 /* Find the parameters. */
1049 wcsdistortion_intermidate_sipeq( cd, cd_inv, a_coeff, b_coeff, pv1, pv2);
1050
1051 /*For a check.
1052 printf("%.10lf\t%.10lf\t%.10lf\t%.10lf\n", cd_inv[0][0], cd_inv[0][1], \
1053 cd_inv[1][0], cd_inv[1][1]);
1054 printf("%.10lf\n", determinant);
1055 */
1056 }
1057
1058
1059
1060
1061
1062 /* Calculate the SIP coefficients from CD matrix parameters and
1063 PV coefficients. */
1064 static double
wcsdistortion_calcsip(size_t axis,size_t m,size_t n,double tpvu[8][8],double tpvv[8][8])1065 wcsdistortion_calcsip(size_t axis, size_t m, size_t n, double tpvu[8][8],
1066 double tpvv[8][8])
1067 {
1068 double sip_coeff;
1069 if( axis == 1) sip_coeff=tpvu[m][n];
1070 else if(axis == 2) sip_coeff=tpvv[m][n];
1071 else error(EXIT_FAILURE, 0, "%s: axis %zu does not exists! ",
1072 __func__, axis);
1073
1074 if( (axis == 1) && (m == 1) && (n == 0) ) sip_coeff = sip_coeff - 1.0;
1075 else if( (axis == 2) && (m == 0) && (n == 1) ) sip_coeff = sip_coeff - 1.0;
1076
1077 return sip_coeff;
1078 }
1079
1080
1081
1082
1083
1084 /* To calculate reverse sip coefficients, we make a grid and polpulate
1085 it with forward coefficients. Then we fit a polynomial through these
1086 points and finally extract coefficeint from these polynomial which are
1087 the reverse sip cefficents. */
1088 static void
wcsdistortion_fitreverse(double * u,double * v,size_t a_order,size_t b_order,size_t naxis1,size_t naxis2,double a_coeff[5][5],double b_coeff[5][5],double ap_coeff[5][5],double bp_coeff[5][5])1089 wcsdistortion_fitreverse(double *u, double *v, size_t a_order, size_t b_order,
1090 size_t naxis1, size_t naxis2, double a_coeff[5][5],
1091 double b_coeff[5][5], double ap_coeff[5][5],
1092 double bp_coeff[5][5])
1093 {
1094 double chisq_ap, chisq_bp;
1095 double *udiff=NULL, *vdiff=NULL;
1096 double *uprime=NULL, *vprime=NULL;
1097 double **udict=NULL, **vdict=NULL;
1098 size_t tsize=(naxis1/4)*(naxis2/4);
1099 double **updict=NULL, **vpdict=NULL;
1100 gsl_vector *y_ap, *y_bp, *c_ap, *c_bp;
1101 size_t ap_order=a_order, bp_order=b_order;
1102 gsl_matrix *X_ap, *X_bp, *cov_ap, *cov_bp;
1103 size_t i=0, j=0, k=0, p_ap=0, p_bp=0, ij=0;
1104 gsl_multifit_linear_workspace *work_ap, *work_bp;
1105
1106 /* Storage structures used:
1107
1108 updict, vpdict - dictionaries (key-value pairs) with key being an
1109 integer starting from 0 to max(aporder,
1110 bporder)(inclusive) and its values are the uprime,
1111 vprime raised to the powers of corresponding keys for
1112 each key.
1113
1114 udict, vdict - dictionaries (key-value pairs) with key being an
1115 integer starting from 0 to max(aorder,
1116 border)(inclusive) and its values are the u, v raised
1117 to the powers of corresponding keys for each key.
1118
1119 u, v - The 1d representation of 2d grid of all points
1120 strating from -CRPIXi to NAXISi - CRPIXi (with a
1121 stride of 4). CRPIXi is subtracted to bring pixels
1122 in world coordinate system (wcs).
1123
1124 uprime, vprime - The grid (represented internally as a 1d array)
1125 with forward coefficients fitted on them
1126 using the relevant udict, vdict values.
1127
1128 udiff, vdiff - 1d array with the values of uprime, vprim subtracted
1129 from u, v arrays.
1130
1131
1132 In matrix equation AX=B,
1133 For axis 1 - A = transpose of a matrix of udict*vdict
1134 B = udiff
1135 For axis 2 - A = transpose of a matrix of updict*vpdict
1136 B = vdiff
1137 */
1138
1139 /* Allocate updict and vpdict and initialize them. */
1140 updict=malloc((ap_order+1)*sizeof(*updict));
1141 vpdict=malloc((bp_order+1)*sizeof(*vpdict));
1142 for(i=0; i<=wcsdistortion_max(ap_order, bp_order); ++i)
1143 {
1144 updict[i]=malloc(tsize*sizeof(updict));
1145 vpdict[i]=malloc(tsize*sizeof(vpdict));
1146 }
1147
1148 /* Allocate and initialize uprime and vprime. */
1149 uprime=malloc(tsize*sizeof(*uprime));
1150 vprime=malloc(tsize*sizeof(*vprime));
1151 for(i=0; i<tsize; ++i)
1152 {
1153 uprime[i]=u[i];
1154 vprime[i]=v[i];
1155 }
1156
1157 /* Allocate and initialize udict and vdict.*/
1158 udict=malloc((a_order+1)*sizeof(*udict));
1159 vdict=malloc((b_order+1)*sizeof(*vdict));
1160 for(i=0; i<=wcsdistortion_max(a_order, b_order); ++i)
1161 {
1162 udict[i]=malloc(tsize*sizeof(udict));
1163 vdict[i]=malloc(tsize*sizeof(vdict));
1164 }
1165 for(i=0; i<tsize; ++i)
1166 {
1167 udict[0][i]=1;
1168 vdict[0][i]=1;
1169 }
1170
1171
1172 /* Fill the values from the in the dicts. The rows of the
1173 dicts act as a key to achieve a key-value functionality. */
1174 for(i=1; i<=a_order; ++i)
1175 for(j=0; j<tsize; ++j)
1176 udict[i][j]=udict[i-1][j]*u[j];
1177
1178 for(i=1; i<=b_order; ++i)
1179 for(j=0; j<tsize; ++j)
1180 vdict[i][j]=vdict[i-1][j]*v[j];
1181
1182
1183 /* Populating forward coefficients on a grid. */
1184 for(i=0; i<=a_order; ++i)
1185 for(j=0; j<=a_order-i; ++j)
1186 {
1187 for(k=0; k<tsize; ++k)
1188 uprime[k]+=a_coeff[i][j]*udict[i][k]*vdict[j][k];
1189
1190 /* The number of parameters for AP_* coefficients. */
1191 ++p_ap;
1192 }
1193 for(i=0; i<=b_order; ++i)
1194 for(j=0; j<=b_order-i; ++j)
1195 {
1196 for(k=0; k<tsize; ++k)
1197 vprime[k]+=b_coeff[i][j]*udict[i][k]*vdict[j][k];
1198
1199 /* The number of parameters for BP_* coefficients. */
1200 ++p_bp;
1201 }
1202
1203 /*For a check.
1204 for(i=0; i<=a_order; ++i)
1205 for(j=0; j<tsize; ++j)
1206 {
1207 printf("udict[%ld][%ld] = %.8lf\n", i, j, uprime[j]);
1208 printf("u%ld = %.10E\n", i, u[i]);
1209 }
1210 */
1211
1212 /* Now we have a grid populated with forward coeffiecients. Now we fit a
1213 reverse polynomial through points using multiparameter linear least
1214 square fittings.*/
1215
1216 /* Initialize dicts. */
1217 for(i=0; i<tsize; ++i)
1218 {
1219 updict[0][i]=1;
1220 vpdict[0][i]=1;
1221 }
1222
1223 /* Fill the values from the in the dicts. The rows of the
1224 dicts act as a key to achieve a key-value functionality. */
1225 for(i=1; i<=ap_order; ++i)
1226 for(j=0; j<tsize; ++j)
1227 updict[i][j]=updict[i-1][j]*uprime[j];
1228
1229 for(i=1; i<=bp_order; ++i)
1230 for(j=0; j<tsize; ++j)
1231 vpdict[i][j]=vpdict[i-1][j]*vprime[j];
1232
1233 /* Allocate memory for Multi-parameter Linear Regressions. */
1234 X_ap = gsl_matrix_alloc (tsize, p_ap);
1235 X_bp = gsl_matrix_alloc (tsize, p_bp);
1236
1237 y_ap = gsl_vector_alloc (tsize);
1238 y_bp = gsl_vector_alloc (tsize);
1239
1240 c_ap = gsl_vector_alloc (p_ap);
1241 c_bp = gsl_vector_alloc (p_bp);
1242
1243 cov_ap = gsl_matrix_alloc (p_ap, p_ap);
1244 cov_bp = gsl_matrix_alloc (p_bp, p_bp);
1245
1246 /* Allocate and initialize udiff and vdiff. */
1247 udiff=malloc(tsize*sizeof(*udiff));
1248 vdiff=malloc(tsize*sizeof(*vdiff));
1249 for(i=0; i<tsize; ++i)
1250 {
1251 udiff[i]=u[i]-uprime[i];
1252 vdiff[i]=v[i]-vprime[i];
1253
1254 gsl_vector_set (y_ap, i, udiff[i]);
1255 gsl_vector_set (y_bp, i, vdiff[i]);
1256 }
1257
1258 /* Filling up he X_ap matrix in column-wise order. */
1259 for(i=0; i<=ap_order; ++i)
1260 for(j=0; j<=ap_order-i; ++j)
1261 {
1262 for(k=0; k<tsize; ++k)
1263 {
1264 gsl_matrix_set(X_ap, k, ij, updict[i][k]*vpdict[j][k]);
1265
1266 /*For a check.
1267 printf("x_ap[%ld] = %.8lf\n", ij, updict[i][k]*vpdict[j][k]);
1268 */
1269 }
1270 ++ij;
1271 }
1272
1273 ij=0;
1274 /* Filling up he X_bp matrix in column-wise order. */
1275 for(i=0; i<=bp_order; ++i)
1276 for(j=0; j<=bp_order-i; ++j)
1277 {
1278 for(k=0; k<tsize; ++k)
1279 {
1280 gsl_matrix_set(X_bp, k, ij, updict[i][k]*vpdict[j][k]);
1281
1282 /*For a check.
1283 printf("x_bp[%ld] = %.8lf\n", ij, updict[i][k]*vpdict[j][k]);
1284 */
1285 }
1286 ++ij;
1287 }
1288
1289 /* Initialize and do the linear least square fitting. */
1290 work_ap = gsl_multifit_linear_alloc (tsize, p_ap);
1291 work_bp = gsl_multifit_linear_alloc (tsize, p_bp);
1292 gsl_multifit_linear(X_ap, y_ap, c_ap, cov_ap, &chisq_ap, work_ap);
1293 gsl_multifit_linear(X_bp, y_bp, c_bp, cov_bp, &chisq_bp, work_bp);
1294
1295 /* Extract the reverse coefficients. */
1296 p_ap=0;
1297 for(i=0; i<=ap_order; ++i)
1298 for(j=0; j<=ap_order-i; ++j)
1299 {
1300 ap_coeff[i][j]=gsl_vector_get(c_ap, p_ap);
1301
1302 /*For a check.
1303 printf("AP_%ld_%ld = %.8E\n", i, j, ap_coeff[i][j]);
1304 */
1305
1306 ++p_ap;
1307 }
1308 p_bp=0;
1309 for(i=0; i<=bp_order; ++i)
1310 for(j=0; j<=bp_order-i; ++j)
1311 {
1312 bp_coeff[i][j]=gsl_vector_get(c_bp, p_bp);
1313
1314 /*For a check.
1315 printf("BP_%ld_%ld = %.8E\n", i, j, bp_coeff[i][j]);
1316 */
1317
1318 ++p_bp;
1319 }
1320
1321 /*For a check.
1322 for(i=0; i<p_ap; ++i)
1323 for(j=0; j<tsize; ++j)
1324 printf("X[%ld][%ld] = %.8lf\n", i, j, gsl_matrix_get(X_ap, j, i));
1325 */
1326
1327 /* Free the memory allocations.*/
1328 gsl_multifit_linear_free(work_bp);
1329 gsl_multifit_linear_free(work_ap);
1330 gsl_matrix_free(cov_bp);
1331 gsl_matrix_free(cov_ap);
1332 gsl_vector_free(c_bp);
1333 gsl_vector_free(c_ap);
1334 gsl_vector_free(y_bp);
1335 gsl_vector_free(y_ap);
1336 gsl_matrix_free(X_bp);
1337 gsl_matrix_free(X_ap);
1338 free(vdiff);
1339 free(udiff);
1340 free(vprime);
1341 free(uprime);
1342
1343 for(i=0; i<ap_order; ++i) { free(vpdict[i]); free(updict[i]); }
1344 for(i=0; i<a_order; ++i) { free(vdict[i]); free(udict[i]); }
1345 free(vpdict);
1346 free(updict);
1347 free(vdict);
1348 free(udict);
1349 }
1350
1351
1352
1353
1354
1355
1356 /* Calculate the reverse sip coefficients. */
1357 static void
wcsdistortion_get_revkeyvalues(struct wcsprm * wcs,size_t * fitsize,double ap_coeff[5][5],double bp_coeff[5][5])1358 wcsdistortion_get_revkeyvalues(struct wcsprm *wcs, size_t *fitsize,
1359 double ap_coeff[5][5], double bp_coeff[5][5])
1360 {
1361 size_t tsize;
1362 size_t i, j, k;
1363 double *u=NULL, *v=NULL;
1364 size_t a_order=0, b_order=0;
1365 double a_coeff[5][5], b_coeff[5][5];
1366 size_t naxis1=fitsize[1], naxis2=fitsize[0];
1367 double crpix1=wcs->crpix[0], crpix2=wcs->crpix[1];
1368
1369 /* Initialise the 2d matrices. */
1370 tsize=(naxis1/4)*(naxis2/4);
1371 for(i=0;i<5;++i) for(j=0;j<5;++j) {a_coeff[i][j]=0; b_coeff[i][j]=0;}
1372
1373 /* Allocate the size of u,v arrays. */
1374 u=malloc(tsize*sizeof(*u));
1375 v=malloc(tsize*sizeof(*v));
1376
1377 if(u==NULL)
1378 error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `u'",
1379 __func__, tsize*sizeof(*u));
1380 if(v==NULL)
1381 error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `v'",
1382 __func__, tsize*sizeof(*v));
1383
1384
1385 /* Make the grid and bring it's origin to the world's origin*/
1386 k=0; for(i=0; i<naxis2; i+=4) for(j=0; j<naxis1; j+=4) u[k++]=j-crpix1;
1387 k=0; for(i=0; i<naxis2; i+=4) for(j=0; j<naxis1; j+=4) v[k++]=i-crpix2;
1388 /*For a check.
1389 for(i=0; i<tsize; ++i)
1390 printf("u%ld = %.10E\n", i, u[i]);
1391 */
1392
1393 wcsdistortion_get_sipcoeff(wcs, &a_order, &b_order, a_coeff, b_coeff);
1394
1395 wcsdistortion_fitreverse(u, v, a_order, b_order, naxis1, naxis2,
1396 a_coeff, b_coeff, ap_coeff, bp_coeff);
1397
1398 /* Free the memory allocations. */
1399 free(v);
1400 free(u);
1401 }
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422 /**************************************************************/
1423 /********** Writing utilities ************/
1424 /**************************************************************/
1425 /* Make the sip key-cards and add them to fullheader.
1426
1427 Return:
1428 char *fullheader - string of keycards.*/
1429 static char *
wcsdistortion_add_sipkeywords(struct wcsprm * wcs,size_t * fitsize,double tpvu[8][8],double tpvv[8][8],int add_reverse,int * nkeys)1430 wcsdistortion_add_sipkeywords(struct wcsprm *wcs, size_t *fitsize,
1431 double tpvu[8][8], double tpvv[8][8],
1432 int add_reverse, int *nkeys)
1433 {
1434 double val=0;
1435 uint8_t i, j, k=0;
1436 int size = wcs->naxis;
1437 size_t a_order=0, b_order=0;
1438 size_t ap_order=0, bp_order=0;
1439 size_t m, n, num=0, numkey=200;
1440 double ap_coeff[5][5], bp_coeff[5][5];
1441
1442 /* The literal strings are a little longer than necessary because GCC
1443 will complain/warn that size_t can be very long. Although it will
1444 never happen that we need to write 10-decimal numbers in these, but to
1445 have a clean build, we'll just set the literal space to be large
1446 enough to avoid the warnings. */
1447 char *fullheader, fmt[50], sipkey[50], keyaxis[50], pcaxis[50];
1448
1449 /* Initialise the 2d matrices. */
1450 *nkeys = 0;
1451 for(i=0;i<5;++i) for(j=0;j<5;++j) {ap_coeff[i][j]=0; bp_coeff[i][j]=0;}
1452
1453 /* The format for each card. */
1454 sprintf(fmt, "%%-8s= %%20.12E%%50s");
1455
1456 /* Allcate memory for cards. */
1457 fullheader = malloc(numkey*80);
1458 if(fullheader==NULL)
1459 error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for `fullheader'",
1460 __func__, sizeof *fullheader);
1461
1462 /* Add other necessary cards. */
1463 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20d%50s", "WCSAXES",
1464 wcs->naxis, "");
1465
1466 for(i=1; i<=size; ++i)
1467 {
1468 sprintf(keyaxis, "CRPIX%d", i);
1469 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.8lf%50s", keyaxis,
1470 wcs->crpix[i-1], "");
1471 }
1472
1473 for(i=1; i<=size; ++i)
1474 for(j=1; j<=size; ++j)
1475 {
1476 sprintf(pcaxis, "PC%d_%d", i, j);
1477 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", pcaxis,
1478 wcs->pc[k++], "");
1479 }
1480
1481 for(i=1; i<=size; ++i)
1482 {
1483 sprintf(keyaxis, "CDELT%d", i);
1484 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", keyaxis,
1485 wcs->cdelt[i-1], "");
1486 }
1487
1488 for(i=1; i<=size; ++i)
1489 {
1490 sprintf(keyaxis, "CUNIT%d", i);
1491 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", keyaxis,
1492 wcs->cunit[i-1]);
1493 }
1494
1495 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", "CTYPE1",
1496 "'RA---TAN-SIP'");
1497 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", "CTYPE2",
1498 "'DEC--TAN-SIP'");
1499
1500 for(i=1; i<=size; ++i)
1501 {
1502 sprintf(keyaxis, "CRVAL%d", i);
1503 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.10lf%50s", keyaxis,
1504 wcs->crval[i-1], "");
1505 }
1506
1507 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", "LONPOLE",
1508 wcs->lonpole, "");
1509 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", "LATPOLE",
1510 wcs->latpole, "");
1511
1512 #if GAL_CONFIG_HAVE_WCSLIB_MJDREF == 1
1513 for(i=1; i<=size; ++i)
1514 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.1lf%50s", "MJDREFI",
1515 wcs->mjdref[i-1], "");
1516 #endif
1517
1518 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", "RADESYS",
1519 wcs->radesys);
1520
1521 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.1lf%50s", "EQUINOX",
1522 wcs->equinox, "");
1523
1524
1525 for(m=0; m<=4; ++m)
1526 for(n=0; n<=4; ++n)
1527 {
1528 /*For axis = 1*/
1529 val=wcsdistortion_calcsip(1, m, n, tpvu, tpvv);
1530 if(val != 0)
1531 {
1532 /* Make keywords */
1533 sprintf(sipkey, "A_%zu_%zu", m, n);
1534 sprintf(fullheader+(FLEN_CARD-1)*num++, fmt, sipkey, val, "");
1535 a_order=wcsdistortion_max(a_order, wcsdistortion_max(m,n));
1536 }
1537
1538 /*For axis = 2*/
1539 val=wcsdistortion_calcsip(2, m, n, tpvu, tpvv);
1540 if(val != 0)
1541 {
1542 /* Make keywords */
1543 sprintf(sipkey, "B_%zu_%zu", m, n);
1544 sprintf(fullheader+(FLEN_CARD-1)*num++, fmt, sipkey, val, "");
1545 b_order=wcsdistortion_max(b_order, wcsdistortion_max(m,n));
1546 }
1547
1548 }
1549
1550 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20zu%50s", "A_ORDER",
1551 a_order, "");
1552 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20zu%50s", "B_ORDER",
1553 b_order, "");
1554
1555 /* If reverse coefficients are required. */
1556 if( add_reverse )
1557 {
1558 ap_order=a_order;
1559 bp_order=b_order;
1560
1561 wcsdistortion_get_revkeyvalues(wcs, fitsize, ap_coeff, bp_coeff);
1562
1563 for(m=0; m<=ap_order; ++m)
1564 for(n=0; n<=ap_order-m; ++n)
1565 {
1566 /*For axis = 1*/
1567 val=ap_coeff[m][n];
1568 if(val != 0)
1569 {
1570 /* Make keywords */
1571 sprintf(sipkey, "AP_%zu_%zu", m, n);
1572 sprintf(fullheader+(FLEN_CARD-1)*num++, fmt, sipkey,
1573 val, "");
1574 }
1575
1576 /*For axis = 2*/
1577 val=bp_coeff[m][n];
1578 if(val != 0)
1579 {
1580 /* Make keywords */
1581 sprintf(sipkey, "BP_%zu_%zu", m, n);
1582 sprintf(fullheader+(FLEN_CARD-1)*num++, fmt, sipkey,
1583 val, "");
1584 }
1585 }
1586
1587 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20zu%50s", "AP_ORDER",
1588 ap_order, "");
1589 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20zu%50s", "BP_ORDER",
1590 bp_order, "");
1591
1592 }
1593
1594 /*For a check.
1595 printf("%s\n", fullheader);
1596 */
1597 *nkeys = num;
1598 return fullheader;
1599 }
1600
1601
1602
1603
1604
1605
1606 /* Make the pv key-cards and add them to fullheader.
1607
1608 Return:
1609 char *fullheader - string of keycards. */
1610 static char *
wcsdistortion_add_pvkeywords(struct wcsprm * wcs,double * pv1,double * pv2,int * nkeys)1611 wcsdistortion_add_pvkeywords(struct wcsprm *wcs, double *pv1,
1612 double *pv2, int *nkeys)
1613 {
1614
1615 double val=0;
1616 uint8_t i, j, k=0;
1617 int size = wcs->naxis;
1618 size_t m, n, num=0, numkey=100;
1619
1620 /* The literal strings are a little longer than necessary because GCC
1621 will complain/warn that size_t can be very long. Although it will
1622 never happen that we need to write 10-decimal numbers in these, but to
1623 have a clean build, we'll just set the literal space to be large
1624 enough to avoid the warnings. */
1625 char *fullheader, fmt[50], pvkey[50], keyaxis[50], pcaxis[50];
1626
1627 /* Initialize values. */
1628 *nkeys = 0;
1629
1630 /* The format for each card. */
1631 sprintf(fmt, "%%-8s= %%20.12E%%50s");
1632
1633 /* Allcate memory for cards. */
1634 fullheader = malloc(numkey*80);
1635 if(fullheader==NULL)
1636 error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for `fullheader'",
1637 __func__, sizeof *fullheader);
1638
1639 /* Add other necessary cards. */
1640 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20d%50s", "WCSAXES",
1641 wcs->naxis, "");
1642
1643 for(i=1; i<=size; ++i)
1644 {
1645 sprintf(keyaxis, "CRPIX%d", i);
1646 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.8lf%50s", keyaxis,
1647 wcs->crpix[i-1], "");
1648 }
1649
1650 for(i=1; i<=size; ++i)
1651 for(j=1; j<=size; ++j)
1652 {
1653 sprintf(pcaxis, "PC%d_%d", i, j);
1654 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", pcaxis,
1655 wcs->pc[k++], "");
1656 }
1657
1658 for(i=1; i<=size; ++i)
1659 {
1660 sprintf(keyaxis, "CDELT%d", i);
1661 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", keyaxis,
1662 wcs->cdelt[i-1], "");
1663 }
1664
1665 for(i=1; i<=size; ++i)
1666 {
1667 sprintf(keyaxis, "CUNIT%d", i);
1668 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", keyaxis,
1669 wcs->cunit[i-1]);
1670 }
1671
1672 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", "CTYPE1",
1673 "'RA---TPV'");
1674 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", "CTYPE2",
1675 "'DEC--TPV'");
1676
1677 for(i=1; i<=size; ++i)
1678 {
1679 sprintf(keyaxis, "CRVAL%d", i);
1680 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.10lf%50s", keyaxis,
1681 wcs->crval[i-1], "");
1682 }
1683
1684 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", "LONPOLE",
1685 wcs->lonpole, "");
1686 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.17lf%50s", "LATPOLE",
1687 wcs->latpole, "");
1688
1689 #if GAL_CONFIG_HAVE_WCSLIB_MJDREF == 1
1690 for(i=1; i<=size; ++i)
1691 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.1lf%50s", "MJDREFI",
1692 wcs->mjdref[i-1], "");
1693 #endif
1694
1695 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %-70s", "RADESYS",
1696 wcs->radesys);
1697
1698 sprintf(fullheader+(FLEN_CARD-1)*num++, "%-8s= %20.1lf%50s", "EQUINOX",
1699 wcs->equinox, "");
1700
1701 for(m=1; m<=2; ++m)
1702 for(n=0; n<=16; ++n)
1703 {
1704 /*For axis = 1*/
1705 if(m == 1)
1706 {
1707 val=pv1[n];
1708 if(val != 0)
1709 {
1710 /* Make keywords */
1711 sprintf(pvkey, "PV%zu_%zu", m, n);
1712 sprintf(fullheader+(FLEN_CARD-1)*num++, fmt, pvkey,
1713 val, "");
1714 }
1715 }
1716
1717 /*For axis = 2*/
1718 if(m == 2)
1719 {
1720 val=pv2[n];
1721 if(val != 0)
1722 {
1723 /* Make keywords */
1724 sprintf(pvkey, "PV%zu_%zu", m, n);
1725 sprintf(fullheader+(FLEN_CARD-1)*num++, fmt, pvkey,
1726 val, "");
1727 }
1728 }
1729
1730 }
1731
1732 /*For a check.
1733 printf("%s\n", fullheader);
1734 */
1735 *nkeys = num;
1736 return fullheader;
1737
1738 }
1739
1740
1741
1742
1743
1744 /* Set the internal structure and do sanity checks. */
1745 static void
wcsdistortion_set_internalstruct(struct wcsprm * wcs,char * fullheader,size_t fulllen,int status,int nwcs,int sumcheck)1746 wcsdistortion_set_internalstruct(struct wcsprm *wcs, char *fullheader,
1747 size_t fulllen, int status, int nwcs,
1748 int sumcheck)
1749 {
1750 size_t i;
1751
1752 if( wcs == NULL )
1753 {
1754
1755 fprintf(stderr, "\n##################\n"
1756 "WCSLIB Warning: wcspih ERROR %d: %s.\n"
1757 "##################\n",
1758 status, wcs_errmsg[status]);
1759 wcs=NULL; nwcs=0;
1760 }
1761
1762 if(wcs)
1763 {
1764 /* It may happen that the WCS-related keyword values are stored as
1765 strings (they have single-quotes around them). In this case,
1766 WCSLIB will read the CRPIX and CRVAL values as zero. When this
1767 happens do a small check and abort, while informing the user about
1768 the problem. */
1769 sumcheck=0;
1770 for(i=0;i<wcs->naxis;++i)
1771 {sumcheck += (wcs->crval[i]==0.0f) + (wcs->crpix[i]==0.0f);}
1772 if(sumcheck==wcs->naxis*2)
1773 {
1774 /* We only care about the first set of characters in each
1775 80-character row, so we don't need to parse the last few
1776 characters anyway. */
1777 fulllen=strlen(fullheader)-12;
1778 for(i=0;i<fulllen;++i)
1779 if( strncmp(fullheader+i, "CRVAL1 = '", 11) == 0 )
1780 fprintf(stderr, "WARNING: WCS Keyword values are not "
1781 "numbers.\n\n"
1782 "WARNING: The values to the WCS-related keywords are "
1783 "enclosed in single-quotes. In the FITS standard "
1784 "this is how string values are stored, therefore "
1785 "WCSLIB is unable to read them AND WILL PUT ZERO IN "
1786 "THEIR PLACE (creating a wrong WCS in the output). "
1787 "Please update the respective keywords of the input "
1788 "to be numbers (see next line).\n\n"
1789 "WARNING: You can do this with Gnuastro's `astfits' "
1790 "program and the `--update' option. The minimal WCS "
1791 "keywords that need a numerical value are: `CRVAL1', "
1792 "`CRVAL2', `CRPIX1', `CRPIX2', `EQUINOX' and "
1793 "`CD%%_%%' (or `PC%%_%%', where the %% are integers), "
1794 "please see the FITS standard, and inspect your FITS "
1795 "file to identify the full set of keywords that you "
1796 "need correct (for example PV%%_%% keywords).\n\n");
1797 }
1798
1799 /* CTYPE is a mandatory WCS keyword, so if it hasn't been given (its
1800 '\0'), then the headers didn't have a WCS structure. However,
1801 WCSLIB still fills in the basic information (for example the
1802 dimensionality of the dataset). */
1803 if(wcs->ctype[0][0]=='\0')
1804 {
1805 wcsfree(wcs);
1806 wcs=NULL;
1807 nwcs=0;
1808 }
1809 else
1810 {
1811 /* For a check.
1812 printf("flag: %d\n", wcs->flag);
1813 printf("naxis: %d\n", wcs->naxis);
1814 printf("crpix: %f, %f\n", wcs->crpix[0], wcs->crpix[1]);
1815 printf("pc: %f, %f, %f, %f\n", wcs->pc[0], wcs->pc[1],
1816 wcs->pc[2], wcs->pc[3]);
1817 printf("cdelt: %f, %f\n", wcs->cdelt[0], wcs->cdelt[1]);
1818 printf("crval: %f, %f\n", wcs->crval[0], wcs->crval[1]);
1819 printf("cunit: %s, %s\n", wcs->cunit[0], wcs->cunit[1]);
1820 printf("ctype: %s, %s\n", wcs->ctype[0], wcs->ctype[1]);
1821 printf("lonpole: %f\n", wcs->lonpole);
1822 printf("latpole: %f\n", wcs->latpole);
1823 */
1824
1825 /* Set the WCS structure. */
1826 status=wcsset(wcs);
1827 if(status)
1828 {
1829 fprintf(stderr, "\n##################\n"
1830 "WCSLIB Warning: wcsset ERROR %d: %s.\n"
1831 "##################\n",
1832 status, wcs_errmsg[status]);
1833 wcsfree(wcs);
1834 wcs=NULL;
1835 nwcs=0;
1836 }
1837 else
1838 /* A correctly useful WCS is present. When no PC matrix
1839 elements were present in the header, the default PC matrix
1840 (a unity matrix) is used. In this case WCSLIB doesn't set
1841 `altlin' (and gives it a value of 0). In Gnuastro, later on,
1842 we might need to know the type of the matrix used, so in
1843 such a case, we will set `altlin' to 1. */
1844 if(wcs->altlin==0) wcs->altlin=1;
1845 }
1846 }
1847
1848 /* For a check.
1849 wcsprt(wcs);
1850 */
1851
1852 /* Clean up. */
1853 status=0;
1854 if (fits_free_memory(fullheader, &status) )
1855 gal_fits_io_error(status, "problem in freeing the memory used to "
1856 "keep all the headers");
1857 }
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878 /************************************************************************/
1879 /*************** High-level functions *********************/
1880 /************************************************************************/
1881 /* After the actual TPV->SIP conversions are done and the keycards are
1882 created using `wcsdistortion_add_sipkeywords`, we then use `wcspih`
1883 to generate and return the new headers having SIP coefficients.
1884
1885 Parameters:
1886 struct wcsprm *inwcs - The wcs parameters of the input fits file.
1887 size_t *fitsize - The size of the array along each dimension.
1888
1889 Return:
1890 struct wcsprm *outwcs - The transformed wcs parameters in the
1891 sip distortion type.*/
1892 struct wcsprm *
gal_wcsdistortion_tpv_to_sip(struct wcsprm * inwcs,size_t * fitsize)1893 gal_wcsdistortion_tpv_to_sip(struct wcsprm *inwcs,
1894 size_t *fitsize)
1895 {
1896 int ctrl=0; /* Don't report why a keyword wasn't used. */
1897 int nreject=0; /* Number of keywords rejected for syntax. */
1898 double cd[2][2];
1899 size_t i=0, j=0;
1900 size_t fulllen=0;
1901 char *fullheader;
1902 int relax=WCSHDR_all; /* Macro: use all informal WCS extensions. */
1903 int nwcs, sumcheck=0;
1904 int nkeys=0, status=0;
1905 struct wcsprm *outwcs=NULL;
1906 double tpvu[8][8], tpvv[8][8];
1907
1908 /* Initialise the 2d matrices. */
1909 for(i=0; i<2; ++i) for(j=0; j<2; ++j) cd[i][j]=0;
1910 for(i=0; i<8; ++i) for(j=0; j<8; ++j) { tpvu[i][j]=0; tpvv[i][j]=0; }
1911
1912 /* Calculate the pv equations and extract sip coefficients from it. */
1913 wcsdistortion_calc_tpveq(inwcs, cd, tpvu, tpvv);
1914
1915 /* Add the sip keywords. */
1916 fullheader=wcsdistortion_add_sipkeywords(inwcs, fitsize, tpvu, tpvv,
1917 1, &nkeys);
1918
1919 /* WCSlib function to parse the FITS headers. */
1920 status=wcspih(fullheader, nkeys, relax, ctrl, &nreject, &nwcs, &outwcs);
1921
1922 /* Set the internal structure. */
1923 wcsdistortion_set_internalstruct(outwcs, fullheader, fulllen, status,
1924 nwcs, sumcheck);
1925
1926 /* Return the output WCS. */
1927 return outwcs;
1928
1929 }
1930
1931
1932
1933
1934
1935 /* After the actual SIP->TPV conversions are done and the keycards are
1936 created using `wcsdistortion_add_pvkeywords`, we then use `wcspih`
1937 to generate and return the new headers having PV coefficients.
1938
1939 Parameters:
1940 struct wcsprm *inwcs - The wcs parameters of the input fits file.
1941
1942 Return:
1943 struct wcsprm *outwcs - The transformed wcs parameters in the
1944 pv distortion type.*/
1945 struct wcsprm *
gal_wcsdistortion_sip_to_tpv(struct wcsprm * inwcs)1946 gal_wcsdistortion_sip_to_tpv(struct wcsprm *inwcs)
1947 {
1948 int ctrl=0; /* Don't report why a keyword wasn't used. */
1949 int nreject=0; /* Number of keywords rejected for syntax. */
1950 double cd[2][2];
1951 size_t i=0, j=0;
1952 size_t fulllen=0;
1953 char *fullheader;
1954 int nwcs, sumcheck=0;
1955 int nkeys=0, status=0;
1956 int relax=WCSHDR_all; /* Macro: use all informal WCS extensions. */
1957 struct wcsprm *outwcs=NULL;
1958 double pv1[17]={0}, pv2[17]={0};
1959
1960 /* Initialise the 2d matrices. */
1961 for(i=0; i<2; ++i) for(j=0; j<2; ++j) cd[i][j]=0;
1962
1963 /* Calculate the sip equations and extract pv coefficients from it. */
1964 wcsdistortion_calc_sipeq(inwcs, cd, pv1, pv2);
1965
1966 /* Add the pv keywords. */
1967 fullheader=wcsdistortion_add_pvkeywords(inwcs, pv1, pv2, &nkeys);
1968
1969
1970 /* WCSlib function to parse the FITS headers. */
1971 status=wcspih(fullheader, nkeys, relax, ctrl, &nreject, &nwcs, &outwcs);
1972
1973 /* Set the internal structure. */
1974 wcsdistortion_set_internalstruct(outwcs, fullheader, fulllen, status,
1975 nwcs, sumcheck);
1976
1977 return outwcs;
1978 }
1979