1 /*********************************************************************
2 Functions to that only use WCSLIB functionality.
3 This is part of GNU Astronomy Utilities (Gnuastro) package.
4
5 Original author:
6 Mohammad Akhlaghi <mohammad@akhlaghi.org>
7 Contributing author(s):
8 Copyright (C) 2015-2021, Free Software Foundation, Inc.
9
10 Gnuastro is free software: you can redistribute it and/or modify it
11 under the terms of the GNU General Public License as published by the
12 Free Software Foundation, either version 3 of the License, or (at your
13 option) any later version.
14
15 Gnuastro is distributed in the hope that it will be useful, but
16 WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
22 **********************************************************************/
23 #include <config.h>
24
25 #include <time.h>
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 #include <assert.h>
33
34 #include <gsl/gsl_linalg.h>
35 #include <wcslib/wcsmath.h>
36 #include <wcslib/wcsprintf.h>
37
38 #include <gnuastro/wcs.h>
39 #include <gnuastro/tile.h>
40 #include <gnuastro/fits.h>
41 #include <gnuastro/pointer.h>
42 #include <gnuastro/dimension.h>
43 #include <gnuastro/statistics.h>
44 #include <gnuastro/permutation.h>
45
46 #include <gnuastro-internal/checkset.h>
47
48 #if GAL_CONFIG_HAVE_WCSLIB_DIS_H
49 #include <wcslib/dis.h>
50 #include <gnuastro-internal/wcsdistortion.h>
51 #endif
52
53
54
55
56
57 /*************************************************************
58 *********** Read WCS ***********
59 *************************************************************/
60 /* It may happen that both the PC+CDELT and CD matrices are present in a
61 file. But in some cases, they may not result in the same rotation
62 matrix. So we need to let the user know about this problem with the FITS
63 file, and as a default behavior, we'll disable the PC matrix (which also
64 needs a CDELT matrix (that may not have been written). */
65 static void
wcs_read_correct_pc_cd(struct wcsprm * wcs)66 wcs_read_correct_pc_cd(struct wcsprm *wcs)
67 {
68 int removepc=0;
69 size_t i, j, naxis=wcs->naxis;
70 double *cdfrompc=gal_pointer_allocate(GAL_TYPE_FLOAT64, naxis*naxis, 0,
71 __func__, "cdfrompc");
72
73 /* A small sanity check. */
74 if(wcs->cdelt==NULL)
75 error(EXIT_FAILURE, 0, "%s: the WCS structure has no 'cdelt' array, "
76 "please contact us at %s to see what the problem is", __func__,
77 PACKAGE_BUGREPORT);
78
79 /* Multiply the PC matrix with the CDELT matrix. */
80 for(i=0;i<naxis;++i)
81 for(j=0;j<naxis;++j)
82 cdfrompc[i*naxis+j] = wcs->cdelt[i] * wcs->pc[i*naxis+j];
83
84 /* Make sure the file's CD matrix is the same as the CD matrix that is
85 derived from the PC+CDELT matrix above. We'll divide the difference by
86 the samller value and if the result is larger than 1e-6, then we'll
87 consider it different. */
88 for(i=0;i<naxis*naxis;++i)
89 if( fabs( wcs->cd[i] - cdfrompc[i] )
90 / ( fabs(wcs->cd[i]) < fabs(cdfrompc[i])
91 ? fabs(wcs->cd[i])
92 : fabs(cdfrompc[i]) )
93 > 1e-5 )
94 { removepc=1; break; }
95
96 /* If the two matrices are different, then print the warning and remove
97 the PC+CDELT matrices to only keep the CD matrix. */
98 if(removepc)
99 {
100 /* Let the user know that there is a problem in the file. */
101 error(EXIT_SUCCESS, 0, "the WCS structure has both the PC matrix "
102 "and CD matrix. However, the two don't match and there is "
103 "no way to know which one was intended by the creator of "
104 "this file. THIS PROGRAM WILL ASSUME THE CD MATRIX AND "
105 "CONTINUE. BUT THIS MAY BE WRONG! To avoid confusion and "
106 "wrong results, its best to only use one of them in your "
107 "FITS file. You can use Gnuastro's 'astfits' program to "
108 "remove any that you want (please run 'info astfits' for "
109 "more). For example if you want to delete the PC matrix "
110 "you can use this command: 'astfits file.fits --delete=PC1_1 "
111 "--delete=PC1_2 --delete=PC2_1 --delete=PC2_2'");
112
113 /* Set the PC matrix to be equal to the CD matrix, and set the CDELTs
114 to 1. */
115 for(i=0;i<naxis;++i) wcs->cdelt[i] = 1.0f;
116 for(i=0;i<naxis*naxis;++i) wcs->pc[i] = wcs->cd[i];
117 }
118
119 /* Clean up. */
120 free(cdfrompc);
121 }
122
123
124
125
126
127 /* For the TPV, TNX and ZPX distortions, WCSLIB can't deal with the CDELT
128 keys properly and its better to use the CD matrix instead, this function
129 will check for this and return 1 if a CD matrix should be used
130 (over-riding the user's desired matrix if necessary). */
131 static int
wcs_use_cd_for_distortion(struct wcsprm * wcs)132 wcs_use_cd_for_distortion(struct wcsprm *wcs)
133 {
134 return ( wcs
135 && wcs->lin.disseq
136 && ( !strcmp( wcs->lin.disseq->dtype[1], "TPV")
137 || !strcmp(wcs->lin.disseq->dtype[1], "TNX")
138 || !strcmp(wcs->lin.disseq->dtype[1], "ZPX") ) );
139 }
140
141
142
143
144
145 /* Read the WCS information from the header. Unfortunately, WCS lib is
146 not thread safe, so it needs a mutex. In case you are not using
147 multiple threads, just pass a NULL pointer as the mutex.
148
149 After you finish with this WCS, you should free the space with:
150
151 status = wcsvfree(&nwcs,&wcs);
152
153 If the WCS structure is not recognized, then this function will
154 return a NULL pointer for the wcsprm structure and a zero for
155 nwcs. It will also report the fact to the user in stderr.
156
157 ===================================
158 WARNING: wcspih IS NOT THREAD SAFE!
159 ===================================
160 Don't call this function within a thread or use a mutex.
161 */
162 struct wcsprm *
gal_wcs_read_fitsptr(fitsfile * fptr,int linearmatrix,size_t hstartwcs,size_t hendwcs,int * nwcs)163 gal_wcs_read_fitsptr(fitsfile *fptr, int linearmatrix, size_t hstartwcs,
164 size_t hendwcs, int *nwcs)
165 {
166 /* Declaratins: */
167 size_t i, fulllen;
168 int nkeys=0, status=0;
169 int sumcheck, nocomments;
170 struct wcsprm *wcs=NULL;
171 char *fullheader, *to, *from;
172 int fixstatus[NWCSFIX]={0};/* For the various wcsfix checks. */
173 int relax = WCSHDR_all; /* Macro: use all informal WCS extensions. */
174 int ctrl = 0; /* Don't report why a keyword wasn't used. */
175 int nreject = 0; /* Number of keywords rejected for syntax. */
176 int fixctrl = 1; /* Correct non-standard units in wcsfix. */
177 void *fixnaxis = NULL; /* For now disable cylfix() with this */
178 /* (because it depends on image size). */
179
180 /* In case the user has asked to limit the HDU keyword cards to use for
181 WCS reading, also count comment/history/empty lines (so the lines here
182 correspond to the output of 'astfits image.fits -h1'). But if no
183 limitation is requested, avoid those lines to make the processing
184 easier for WCSLIB. */
185 nocomments = hendwcs>hstartwcs ? 0 : 1;
186
187 /* CFITSIO function: */
188 if( fits_hdr2str(fptr, nocomments, NULL, 0, &fullheader,
189 &nkeys, &status) )
190 gal_fits_io_error(status, NULL);
191
192 /* Only consider the header keywords in the current range: */
193 if(hendwcs>hstartwcs)
194 {
195 /* Mark the last character in the desired region. */
196 fullheader[hendwcs*(FLEN_CARD-1)]='\0';
197
198 /* For a check:
199 printf("%s\n", fullheader);
200 */
201
202 /* Shift all the characters to the start of the string. */
203 if(hstartwcs) /* hstartwcs!=0 */
204 {
205 to=fullheader;
206 from=&fullheader[hstartwcs*(FLEN_CARD-1)-1];
207 while(*from++!='\0') *to++=*from;
208 }
209
210 nkeys=hendwcs-hstartwcs;
211
212 /* For a check:
213 printf("\n\n\n###############\n\n\n\n\n\n");
214 printf("%s\n", &fullheader[1*(FLEN_CARD-1)]);
215 exit(0);
216 */
217 }
218
219 /* WCSlib function to parse the FITS headers. */
220 status=wcspih(fullheader, nkeys, relax, ctrl, &nreject, nwcs, &wcs);
221 if(status)
222 {
223 fprintf(stderr,"\n"
224 "##################\n"
225 "WCSLIB Warning: wcspih ERROR %d: %s.\n"
226 "##################\n",
227 status, wcs_errmsg[status]);
228 wcs=NULL; *nwcs=0;
229 }
230
231 /* Set the internal structure: */
232 if(wcs)
233 {
234 /* It may happen that the WCS-related keyword values are stored as
235 strings (they have single-quotes around them). In this case,
236 WCSLIB will read the CRPIX and CRVAL values as zero. When this
237 happens do a small check and abort, while informing the user about
238 the problem. */
239 sumcheck=0;
240 for(i=0;i<wcs->naxis;++i)
241 {sumcheck += (wcs->crval[i]==0.0f) + (wcs->crpix[i]==0.0f);}
242 if(sumcheck==wcs->naxis*2)
243 {
244 /* We only care about the first set of characters in each
245 80-character row, so we don't need to parse the last few
246 characters anyway. */
247 fulllen=strlen(fullheader)-12;
248 for(i=0;i<fulllen;++i)
249 if( strncmp(fullheader+i, "CRVAL1 = '", 11) == 0 )
250 fprintf(stderr, "WARNING: WCS Keyword values are not "
251 "numbers.\n\n"
252 "WARNING: The values to the WCS-related keywords are "
253 "enclosed in single-quotes. In the FITS standard "
254 "this is how string values are stored, therefore "
255 "WCSLIB is unable to read them AND WILL PUT ZERO IN "
256 "THEIR PLACE (creating a wrong WCS in the output). "
257 "Please update the respective keywords of the input "
258 "to be numbers (see next line).\n\n"
259 "WARNING: You can do this with Gnuastro's 'astfits' "
260 "program and the '--update' option. The minimal WCS "
261 "keywords that need a numerical value are: 'CRVAL1', "
262 "'CRVAL2', 'CRPIX1', 'CRPIX2', 'EQUINOX' and "
263 "'CD%%_%%' (or 'PC%%_%%', where the %% are integers), "
264 "please see the FITS standard, and inspect your FITS "
265 "file to identify the full set of keywords that you "
266 "need correct (for example PV%%_%% keywords).\n\n");
267 }
268
269 /* CTYPE is a mandatory WCS keyword, so if it hasn't been given (its
270 '\0'), then the headers didn't have a WCS structure. However,
271 WCSLIB still fills in the basic information (for example the
272 dimensionality of the dataset). */
273 if(wcs->ctype[0][0]=='\0')
274 {
275 wcsfree(wcs);
276 wcs=NULL;
277 *nwcs=0;
278 }
279 else
280 {
281 /* For a check (we can't use 'wcsprt(wcs)' because this WCS isn't
282 yet initialized).
283 printf("flag: %d\n", wcs->flag);
284 printf("NAXIS: %d\n", wcs->naxis);
285 printf("CRPIX: ");
286 for(i=0;i<wcs->naxis;++i)
287 { printf("%g, ", wcs->crpix[i]); } printf("\n");
288 printf("PC: ");
289 for(i=0;i<wcs->naxis*wcs->naxis;++i)
290 { printf("%g, ", wcs->pc[i]); } printf("\n");
291 printf("CDELT: ");
292 for(i=0;i<wcs->naxis;++i)
293 { printf("%g, ", wcs->cdelt[i]);} printf("\n");
294 printf("CD: ");
295 for(i=0;i<wcs->naxis*wcs->naxis;++i)
296 { printf("%g, ", wcs->cd[i]); } printf("\n");
297 printf("CRVAL: ");
298 for(i=0;i<wcs->naxis;++i)
299 { printf("%g, ", wcs->crval[i]); } printf("\n");
300 printf("CUNIT: ");
301 for(i=0;i<wcs->naxis;++i)
302 { printf("%s, ", wcs->cunit[i]); } printf("\n");
303 printf("CTYPE: ");
304 for(i=0;i<wcs->naxis;++i)
305 { printf("%s, ", wcs->ctype[i]); } printf("\n");
306 printf("LONPOLE: %f\n", wcs->lonpole);
307 printf("LATPOLE: %f\n", wcs->latpole);
308 */
309
310 /* Some datasets may use 'angstroms' (not case-sensitive) in the
311 third dimension instead of the standard 'angstrom' (note the
312 differing 's'). In this case WCSLIB (atleast until version
313 7.3) will not recognize it. We will therefore manually remove
314 the 's' before feeding the WCS structure to WCSLIB. */
315 if( wcs->naxis==3
316 && strlen(wcs->cunit[2])==9
317 && !strncasecmp(wcs->cunit[2], "angstroms", 9) )
318 wcs->cunit[2][8]='\0';
319
320 /* Fix non-standard WCS features. */
321 if( wcsfix(fixctrl, fixnaxis, wcs, fixstatus) )
322 {
323 if(fixstatus[CDFIX])
324 error(0, 0, "%s: (warning) wcsfix status for cdfix: %d",
325 __func__, fixstatus[CDFIX]);
326 if(fixstatus[DATFIX])
327 error(0, 0, "%s: (warning) wcsfix status for datfix: %d",
328 __func__, fixstatus[DATFIX]);
329 #if GAL_CONFIG_HAVE_WCSLIB_OBSFIX
330 if(fixstatus[OBSFIX])
331 error(0, 0, "%s: (warning) wcsfix status for obsfix: %d",
332 __func__, fixstatus[OBSFIX]);
333 #endif
334 if(fixstatus[UNITFIX])
335 error(0, 0, "%s: (warning) wcsfix status for unitfix: %d",
336 __func__, fixstatus[UNITFIX]);
337 if(fixstatus[SPCFIX])
338 error(0, 0, "%s: (warning) wcsfix status for spcfix: %d",
339 __func__, fixstatus[SPCFIX]);
340 if(fixstatus[CELFIX])
341 error(0, 0, "%s: (warning) wcsfix status for celfix: %d",
342 __func__, fixstatus[CELFIX]);
343 if(fixstatus[CYLFIX])
344 error(0, 0, "%s: (warning) wcsfix status for cylfix: %d",
345 __func__, fixstatus[CYLFIX]);
346 }
347
348 /* Enable wcserr if necessary (can help in debugging situations
349 when the default error message isn't enough). If you confront
350 an error, un-comment the two lines below and put
351 'wcsperr(wcs,NULL);' after the problematic WCSLIB
352 command. This will print a trace-back of the cause.
353 wcserr_enable(1);
354 wcsprintf_set(stderr);
355 */
356
357 /* Set the WCS structure. */
358 status=wcsset(wcs);
359 if(status)
360 {
361 fprintf(stderr, "\n##################\n"
362 "WCSLIB Warning: wcsset ERROR %d: %s.\n"
363 "##################\n",
364 status, wcs_errmsg[status]);
365 wcsfree(wcs);
366 wcs=NULL;
367 *nwcs=0;
368 }
369 /* A correctly useful WCS is present. */
370 else
371 {
372 /* According to WCSLIB in discussing 'altlin': "If none of
373 these bits are set the PCi_ja representation results, i.e.
374 wcsprm::pc and wcsprm::cdelt will be used as given". In
375 effect it will also set the PC matrix to unity. So we can
376 safely set it to '1' here because some parts of Gnuastro
377 will later look into this. */
378 if(wcs->altlin==0) wcs->altlin=1;
379
380 /* If both the PC and CD matrix have been given, the first
381 two bits of 'altlin' will be '1'. We need to make sure
382 they are the same matrix, and let the user know if they
383 aren't. */
384 if( (wcs->altlin & 1) && (wcs->altlin & 2) )
385 wcs_read_correct_pc_cd(wcs);
386 }
387 }
388 }
389
390 /* If the distortion requires a CD matrix, or if the user wants it, do
391 the conversion here, otherwise, make sure the PC matrix is used. */
392 if(wcs_use_cd_for_distortion(wcs) \
393 || linearmatrix==GAL_WCS_LINEAR_MATRIX_CD)
394 gal_wcs_to_cd(wcs);
395 else gal_wcs_decompose_pc_cdelt(wcs);
396
397 /* Clean up and return. */
398 status=0;
399 if (fits_free_memory(fullheader, &status) )
400 gal_fits_io_error(status, "problem in freeing the memory used to "
401 "keep all the headers");
402 return wcs;
403 }
404
405
406
407
408
409 struct wcsprm *
gal_wcs_read(char * filename,char * hdu,int linearmatrix,size_t hstartwcs,size_t hendwcs,int * nwcs)410 gal_wcs_read(char *filename, char *hdu, int linearmatrix,
411 size_t hstartwcs, size_t hendwcs, int *nwcs)
412 {
413 int status=0;
414 fitsfile *fptr;
415 struct wcsprm *wcs;
416
417 /* Make sure we are dealing with a FITS file. */
418 if( gal_fits_file_recognized(filename) == 0 )
419 return NULL;
420
421 /* Check HDU for realistic conditions: */
422 fptr=gal_fits_hdu_open_format(filename, hdu, 0);
423
424 /* Read the WCS information: */
425 wcs=gal_wcs_read_fitsptr(fptr, linearmatrix, hstartwcs,
426 hendwcs, nwcs);
427
428 /* Close the FITS file and return. */
429 fits_close_file(fptr, &status);
430 gal_fits_io_error(status, NULL);
431 return wcs;
432 }
433
434
435
436
437
438 struct wcsprm *
gal_wcs_create(double * crpix,double * crval,double * cdelt,double * pc,char ** cunit,char ** ctype,size_t ndim,int linearmatrix)439 gal_wcs_create(double *crpix, double *crval, double *cdelt,
440 double *pc, char **cunit, char **ctype,
441 size_t ndim, int linearmatrix)
442 {
443 size_t i;
444 int status;
445 struct wcsprm *wcs;
446 double equinox=2000.0f;
447
448 /* Allocate the memory necessary for the wcsprm structure. */
449 errno=0;
450 wcs=malloc(sizeof *wcs);
451 if(wcs==NULL)
452 error(EXIT_FAILURE, errno, "%zu for wcs in preparewcs", sizeof *wcs);
453
454 /* Initialize the structure (allocate all its internal arrays). */
455 wcs->flag=-1;
456 if( (status=wcsini(1, ndim, wcs)) )
457 error(EXIT_FAILURE, 0, "wcsini error %d: %s",
458 status, wcs_errmsg[status]);
459
460 /* Fill in all the important WCS structure parameters. */
461 wcs->altlin = 0x1;
462 wcs->equinox = equinox;
463 for(i=0;i<ndim;++i)
464 {
465 wcs->crpix[i] = crpix[i];
466 wcs->crval[i] = crval[i];
467 wcs->cdelt[i] = cdelt[i];
468 if(cunit[i]) strcpy(wcs->cunit[i], cunit[i]);
469 if(ctype[i]) strcpy(wcs->ctype[i], ctype[i]);
470 }
471 for(i=0;i<ndim*ndim;++i) wcs->pc[i]=pc[i];
472
473 /* Set up the wcs structure with the constants defined above. */
474 status=wcsset(wcs);
475 if(status)
476 error(EXIT_FAILURE, 0, "wcsset error %d: %s", status,
477 wcs_errmsg[status]);
478
479 /* If a CD matrix is desired make it. */
480 if(linearmatrix==GAL_WCS_LINEAR_MATRIX_CD)
481 gal_wcs_to_cd(wcs);
482
483 /* Return the output WCS. */
484 return wcs;
485 }
486
487
488
489
490
491 /* Extract the dimension name from CTYPE. */
492 char *
gal_wcs_dimension_name(struct wcsprm * wcs,size_t dimension)493 gal_wcs_dimension_name(struct wcsprm *wcs, size_t dimension)
494 {
495 size_t i;
496 char *out;
497
498 /* Make sure a WCS pointer actually exists. */
499 if(wcs==NULL) return NULL;
500
501 /* Make sure the requested dimension is not larger than the number of
502 dimensions in the WCS. */
503 if(dimension >= wcs->naxis) return NULL;
504
505 /* Make a copy of the CTYPE value and set the first occurance of '-' to
506 '\0', to avoid the projection type. */
507 gal_checkset_allocate_copy(wcs->ctype[dimension], &out);
508 for(i=0;i<strlen(out);++i) if(out[i]=='-') out[i]='\0';
509
510 /* Return the output array. */
511 return out;
512 }
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532 /*************************************************************
533 *********** Write WCS ***********
534 *************************************************************/
535 char *
gal_wcs_write_wcsstr(struct wcsprm * wcs,int * nkeyrec)536 gal_wcs_write_wcsstr(struct wcsprm *wcs, int *nkeyrec)
537 {
538 char *wcsstr;
539 int status=0;
540
541 /* Finalize the linear transformation matrix. Note that some programs may
542 have worked on the WCS. So even if 'altlin' is already 2, we'll just
543 ensure that the final matrix is CD here. */
544 if(wcs->altlin==2 || wcs_use_cd_for_distortion(wcs)) gal_wcs_to_cd(wcs);
545 else gal_wcs_decompose_pc_cdelt(wcs);
546
547 /* Clean up small errors in the PC matrix and CDELT values. */
548 gal_wcs_clean_errors(wcs);
549
550 /* Convert the WCS information to text. If there was an error, then free
551 any allocated space and put zero on 'nkeyrec'. */
552 status=wcshdo(WCSHDO_safe, wcs, nkeyrec, &wcsstr);
553 if(status)
554 {
555 error(0, 0, "%s: WARNING: WCSLIB error, no WCS in output.\n"
556 "wcshdo ERROR %d: %s", __func__, status, wcs_errmsg[status]);
557 if(wcsstr) free(wcsstr);
558 *nkeyrec=0;
559 wcsstr=NULL;
560 }
561
562 /* Return the string. */
563 return wcsstr;
564 }
565
566
567
568
569
570 void
gal_wcs_write_in_fitsptr(fitsfile * fptr,struct wcsprm * wcs)571 gal_wcs_write_in_fitsptr(fitsfile *fptr, struct wcsprm *wcs)
572 {
573 char *wcsstr;
574 int nkeyrec=0;
575
576 /* Convert the WCS structure into FITS keywords as a string. */
577 wcsstr=gal_wcs_write_wcsstr(wcs, &nkeyrec);
578
579 /* Write the keywords into the FITS file. */
580 gal_fits_key_write_wcsstr(fptr, wcs, wcsstr, nkeyrec);
581 free(wcsstr);
582 }
583
584
585
586
587
588 void
gal_wcs_write(struct wcsprm * wcs,char * filename,char * extname,gal_fits_list_key_t * headers,char * program_string)589 gal_wcs_write(struct wcsprm *wcs, char *filename,
590 char *extname, gal_fits_list_key_t *headers,
591 char *program_string)
592 {
593 int status=0;
594 size_t ndim=0;
595 fitsfile *fptr;
596 long *naxes=NULL;
597
598 /* Small sanity checks */
599 if(wcs==NULL)
600 error(EXIT_FAILURE, 0, "%s: input WCS is NULL", __func__);
601 if( gal_fits_file_recognized(filename)==0 )
602 error(EXIT_FAILURE, 0, "%s: not a FITS suffix", filename);
603
604 /* Open the file for writing */
605 fptr=gal_fits_open_to_write(filename);
606
607 /* Create the FITS file. */
608 fits_create_img(fptr, gal_fits_type_to_bitpix(GAL_TYPE_UINT8),
609 ndim, naxes, &status);
610 gal_fits_io_error(status, NULL);
611
612 /* Remove the two comment lines put by CFITSIO. Note that in some cases,
613 it might not exist. When this happens, the status value will be
614 non-zero. We don't care about this error, so to be safe, we will just
615 reset the status variable after these calls. */
616 fits_delete_key(fptr, "COMMENT", &status);
617 fits_delete_key(fptr, "COMMENT", &status);
618 status=0;
619
620 /* If an extension name was requested, add it. */
621 if(extname)
622 fits_write_key(fptr, TSTRING, "EXTNAME", extname, "", &status);
623
624 /* Write the WCS structure. */
625 gal_wcs_write_in_fitsptr(fptr, wcs);
626
627 /* Write all the headers and the version information. */
628 gal_fits_key_write_version_in_ptr(&headers, program_string, fptr);
629
630 /* Close the FITS file. */
631 fits_close_file(fptr, &status);
632 gal_fits_io_error(status, NULL);
633 }
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654 /*************************************************************
655 *********** Coordinate system ***********
656 *************************************************************/
657 int
gal_wcs_coordsys_from_string(char * coordsys)658 gal_wcs_coordsys_from_string(char *coordsys)
659 {
660 if( !strcmp(coordsys,"eq-j2000") ) return GAL_WCS_COORDSYS_EQJ2000;
661 else if( !strcmp(coordsys,"eq-b1950") ) return GAL_WCS_COORDSYS_EQB1950;
662 else if( !strcmp(coordsys,"ec-j2000") ) return GAL_WCS_COORDSYS_ECJ2000;
663 else if( !strcmp(coordsys,"ec-b1950") ) return GAL_WCS_COORDSYS_ECB1950;
664 else if( !strcmp(coordsys,"galactic") ) return GAL_WCS_COORDSYS_GALACTIC;
665 else if( !strcmp(coordsys,"supergalactic") )
666 return GAL_WCS_COORDSYS_SUPERGALACTIC;
667 else
668 error(EXIT_FAILURE, 0, "WCS coordinate system name '%s' not "
669 "recognized, currently recognized names are 'eq-j2000', "
670 "'eq-b1950', 'galactic' and 'supergalactic'", coordsys);
671
672 /* Control should not reach here. */
673 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
674 "problem. Control should not reach the end of this function",
675 __func__, PACKAGE_BUGREPORT);
676 return GAL_WCS_COORDSYS_INVALID;
677 }
678
679
680
681
682 /* Identify the coordinate system of the WCS. */
683 int
gal_wcs_coordsys_identify(struct wcsprm * wcs)684 gal_wcs_coordsys_identify(struct wcsprm *wcs)
685 {
686 /* Equatorial (we are keeping the dash ('-') to make sure it is a
687 standard). */
688 if ( !strncmp(wcs->ctype[0], "RA---", 5)
689 && !strncmp(wcs->ctype[1], "DEC--", 5) )
690 {
691 if ( !strncmp(wcs->radesys, "FK4", 3) )
692 return GAL_WCS_COORDSYS_EQB1950;
693 else if ( !strncmp(wcs->radesys, "FK5", 3) )
694 return GAL_WCS_COORDSYS_EQJ2000;
695 else
696 error(EXIT_FAILURE, 0, "%s: the '%s' value for 'RADESYS' is "
697 "not yet implemented! Please contact us at %s to "
698 "implement it", __func__, wcs->radesys, PACKAGE_BUGREPORT);
699 }
700
701 /* Ecliptic. */
702 else if ( !strncmp(wcs->ctype[0], "ELON-", 5)
703 && !strncmp(wcs->ctype[1], "ELAT-", 5) )
704 if ( !strncmp(wcs->radesys, "FK4", 3) )
705 return GAL_WCS_COORDSYS_ECB1950;
706 else if ( !strncmp(wcs->radesys, "FK5", 3) )
707 return GAL_WCS_COORDSYS_ECJ2000;
708 else
709 error(EXIT_FAILURE, 0, "%s: the '%s' value for 'RADESYS' is "
710 "not yet implemented! Please contact us at %s to "
711 "implement it", __func__, wcs->radesys, PACKAGE_BUGREPORT);
712
713 /* Galactic. */
714 else if ( !strncmp(wcs->ctype[0], "GLON-", 5)
715 && !strncmp(wcs->ctype[1], "GLAT-", 5) )
716 return GAL_WCS_COORDSYS_GALACTIC;
717
718 /* SuperGalactic. */
719 else if ( !strncmp(wcs->ctype[0], "SLON-", 5)
720 && !strncmp(wcs->ctype[1], "SLAT-", 5) )
721 return GAL_WCS_COORDSYS_SUPERGALACTIC;
722
723 /* Other. */
724 else
725 error(EXIT_FAILURE, 0, "%s: the CTYPE values '%s' and '%s' are "
726 "not yet implemented! Please contact us at %s to "
727 "implement it", __func__, wcs->ctype[0], wcs->ctype[1],
728 PACKAGE_BUGREPORT);
729
730 /* Control should not reach here. */
731 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
732 "problem. Control should not reach the end of this function",
733 __func__, PACKAGE_BUGREPORT);
734 return GAL_WCS_COORDSYS_INVALID;
735 }
736
737
738
739
740
741 /* Set the pole coordinates (current values taken from the WCSLIB
742 manual.
743 lng2p1: pole of input (1) system in output (2) system's logitude.
744 lat2p1: pole of input (1) system in output (2) system's latitude.
745 lng1p2: pole of output (2) system in input (1) system's longitude.
746
747 Values from NED (inspired by WCSLIB manual's example).
748 https://ned.ipac.caltech.edu/coordinate_calculator
749
750 longi (deg) latit (deg) OUTPUT INPUT
751 ----- ----- ------ -----
752 (------------, -----------) B1950 equ. coords. of B1950 equ. pole.
753 (180.31684301, 89.72174782) J2000 equ. coords. of B1950 equ. pole.
754 (90.000000000, 66.55421111) B1950 ecl. coords. of B1950 equ. pole.
755 (90.699521110, 66.56068919) J2000 ecl. coords. of B1950 equ. pole.
756 (123.00000000, 27.40000000) Galactic coords. of B1950 equ. pole.
757 (26.731537070, 15.64407736) Supgalactic coords. of B1950 equ. pole.
758
759 (359.68621044, 89.72178502) B1950 equ. coords. of J2000 equ. pole.
760 (------------, -----------) J2000 equ. coords. of J2000 equ. pole.
761 (89.300755510, 66.55417728) B1950 ecl. coords. of J2000 equ. pole.
762 (90.000000000, 66.56070889) J2000 ecl. coords. of J2000 equ. pole.
763 (122.93200023, 27.12843056) Galactic coords. of J2000 equ. pole.
764 (26.450516650, 15.70886131) Supgalactic coords. of J2000 equ. pole.
765
766 (270.00000000, 66.55421111) B1950 equ. coords. of B1950 ecl. pole.
767 (269.99920697, 66.55421892) J2000 equ. coords. of B1950 ecl. pole.
768 (------------, -----------) B1950 ecl. coords. of B1950 ecl. pole.
769 (267.21656404, 89.99350237) J2000 ecl. coords. of B1950 ecl. pole.
770 (96.376479150, 29.81195400) Galactic coords. of B1950 ecl. pole.
771 (33.378919140, 38.34766498) Supgalactic coords. of B1950 ecl. pole.
772
773 (270.00099211, 66.56069675) B1950 equ. coords. of J2000 ecl. pole.
774 (270.00000000, 66.56070889) J2000 equ. coords. of J2000 ecl. pole.
775 (86.517962160, 89.99350236) B1950 ecl. coords. of J2000 ecl. pole.
776 (------------, -----------) J2000 ecl. coords. of J2000 ecl. pole.
777 (96.383958840, 29.81163604) Galactic coords. of J2000 ecl. pole.
778 (33.376119480, 38.34154959) Supgalactic coords. of J2000 ecl. pole.
779
780 (192.25000000, 27.40000000) B1950 equ. coords. of Galactic pole.
781 (192.85949646, 27.12835323) J2000 equ. coords. of Galactic pole.
782 (179.32094769, 29.81195400) B1950 ecl. coords. of Galactic pole.
783 (180.02317894, 29.81153742) J2000 ecl. coords. of Galactic pole.
784 (------------, -----------) Galactic coords. of Galactic pole.
785 (90.000000000, 6.320000000) Supgalactic coords. of Galactic pole.
786
787 (283.18940711, 15.64407736) B1950 equ. coords. of SupGalactic pole.
788 (283.75420420, 15.70894043) J2000 equ. coords. of SupGalactic pole.
789 (286.26975051, 38.34766498) B1950 ecl. coords. of SupGalactic pole.
790 (286.96654469, 38.34158720) J2000 ecl. coords. of SupGalactic pole.
791 (47.370000000, 6.320000000) Galactic coords. of SupGalactic pole.
792 (------------, -----------) Supgalactic coords. of SupGalactic pole.
793 */
794 static void
wcs_coordsys_insys_pole_in_outsys(int insys,int outsys,double * lng2p1,double * lat2p1,double * lng1p2)795 wcs_coordsys_insys_pole_in_outsys(int insys, int outsys, double *lng2p1,
796 double *lat2p1, double *lng1p2)
797 {
798 switch( insys )
799 {
800 case GAL_WCS_COORDSYS_EQB1950:
801 switch( outsys)
802 {
803 case GAL_WCS_COORDSYS_EQB1950:
804 *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN; return;
805 case GAL_WCS_COORDSYS_EQJ2000:
806 *lng2p1=180.31684301; *lat2p1=89.72174782; *lng1p2=359.68621044; return;
807 case GAL_WCS_COORDSYS_ECB1950:
808 *lng2p1=90.000000000; *lat2p1=66.55421111; *lng1p2=270.00000000; return;
809 case GAL_WCS_COORDSYS_ECJ2000:
810 *lng2p1=90.699521110; *lat2p1=66.56068919; *lng1p2=270.00099211; return;
811 case GAL_WCS_COORDSYS_GALACTIC:
812 *lng2p1=123.00000000; *lat2p1=27.40000000; *lng1p2=192.25000000; return;
813 case GAL_WCS_COORDSYS_SUPERGALACTIC:
814 *lng2p1=26.731537070; *lat2p1=15.64407736; *lng1p2=283.18940711; return;
815 default:
816 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
817 "fix the problem. The code '%d' isn't a recognized WCS "
818 "coordinate system ID for 'outsys' (input EQB1950)", __func__,
819 PACKAGE_BUGREPORT, outsys);
820 }
821 break;
822 case GAL_WCS_COORDSYS_EQJ2000:
823 switch( outsys)
824 {
825 case GAL_WCS_COORDSYS_EQB1950:
826 *lng2p1=359.68621044; *lat2p1=89.72178502; *lng1p2=180.31684301; return;
827 case GAL_WCS_COORDSYS_EQJ2000:
828 *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN; return;
829 case GAL_WCS_COORDSYS_ECB1950:
830 *lng2p1=89.300755510; *lat2p1=66.55417728; *lng1p2=269.99920697; return;
831 case GAL_WCS_COORDSYS_ECJ2000:
832 *lng2p1=90.000000000; *lat2p1=66.56070889; *lng1p2=270.00000000; return;
833 case GAL_WCS_COORDSYS_GALACTIC:
834 *lng2p1=122.93200023; *lat2p1=27.12843056; *lng1p2=192.85949646; return;
835 case GAL_WCS_COORDSYS_SUPERGALACTIC:
836 *lng2p1=26.450516650; *lat2p1=15.70886131; *lng1p2=283.75420420; return;
837 default:
838 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
839 "fix the problem. The code '%d' isn't a recognized WCS "
840 "coordinate system ID for 'outsys' (input EQJ2000)", __func__,
841 PACKAGE_BUGREPORT, outsys);
842 }
843 break;
844 case GAL_WCS_COORDSYS_ECB1950:
845 switch( outsys)
846 {
847 case GAL_WCS_COORDSYS_EQB1950:
848 *lng2p1=270.00000000; *lat2p1=66.55421111; *lng1p2=90.000000000; return;
849 case GAL_WCS_COORDSYS_EQJ2000:
850 *lng2p1=269.99920697; *lat2p1=66.55421892; *lng1p2=89.300755510; return;
851 case GAL_WCS_COORDSYS_ECB1950:
852 *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN; return;
853 case GAL_WCS_COORDSYS_ECJ2000:
854 *lng2p1=267.21656404; *lat2p1=89.99350237; *lng1p2=86.517962160; return;
855 case GAL_WCS_COORDSYS_GALACTIC:
856 *lng2p1=96.383958840; *lat2p1=29.81163604; *lng1p2=179.32094769; return;
857 case GAL_WCS_COORDSYS_SUPERGALACTIC:
858 *lng2p1=33.378919140; *lat2p1=38.34766498; *lng1p2=286.26975051; return;
859 default:
860 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
861 "fix the problem. The code '%d' isn't a recognized WCS "
862 "coordinate system ID for 'outsys' (input ECB1950)", __func__,
863 PACKAGE_BUGREPORT, outsys);
864 }
865 break;
866 case GAL_WCS_COORDSYS_ECJ2000:
867 switch( outsys)
868 {
869 case GAL_WCS_COORDSYS_EQB1950:
870 *lng2p1=270.00099211; *lat2p1=66.56069675; *lng1p2=90.699521110; return;
871 case GAL_WCS_COORDSYS_EQJ2000:
872 *lng2p1=270.00000000; *lat2p1=66.56070889; *lng1p2=90.000000000; return;
873 case GAL_WCS_COORDSYS_ECB1950:
874 *lng2p1=86.517962160; *lat2p1=89.99350236; *lng1p2=267.21656404; return;
875 case GAL_WCS_COORDSYS_ECJ2000:
876 *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN; return;
877 case GAL_WCS_COORDSYS_GALACTIC:
878 *lng2p1=96.383958840; *lat2p1=29.81163604; *lng1p2=180.02317894; return;
879 case GAL_WCS_COORDSYS_SUPERGALACTIC:
880 *lng2p1=33.376119480; *lat2p1=38.34154959; *lng1p2=286.96654469; return;
881 default:
882 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
883 "fix the problem. The code '%d' isn't a recognized WCS "
884 "coordinate system ID for 'outsys' (input ECJ2000)", __func__,
885 PACKAGE_BUGREPORT, outsys);
886 }
887 break;
888 case GAL_WCS_COORDSYS_GALACTIC:
889 switch( outsys)
890 {
891 case GAL_WCS_COORDSYS_EQB1950:
892 *lng2p1=192.25000000; *lat2p1=27.40000000; *lng1p2=123.00000000; return;
893 case GAL_WCS_COORDSYS_EQJ2000:
894 *lng2p1=192.85949646; *lat2p1=27.12835323; *lng1p2=122.93200023; return;
895 case GAL_WCS_COORDSYS_ECB1950:
896 *lng2p1=179.32094769; *lat2p1=29.81195400; *lng1p2=96.376479150; return;
897 case GAL_WCS_COORDSYS_ECJ2000:
898 *lng2p1=180.02317894; *lat2p1=29.81153742; *lng1p2=96.383958840; return;
899 case GAL_WCS_COORDSYS_GALACTIC:
900 *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN; return;
901 case GAL_WCS_COORDSYS_SUPERGALACTIC:
902 *lng2p1=90.000000000; *lat2p1=6.320000000; *lng1p2=47.370000000; return;
903 default:
904 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
905 "fix the problem. The code '%d' isn't a recognized WCS "
906 "coordinate system ID for 'outsys' (input GALACTIC)", __func__,
907 PACKAGE_BUGREPORT, outsys);
908 }
909 break;
910 case GAL_WCS_COORDSYS_SUPERGALACTIC:
911 switch( outsys)
912 {
913 case GAL_WCS_COORDSYS_EQB1950:
914 *lng2p1=283.18940711; *lat2p1=15.64407736; *lng1p2=26.731537070; return;
915 case GAL_WCS_COORDSYS_EQJ2000:
916 *lng2p1=283.75420420; *lat2p1=15.70894043; *lng1p2=26.450516650; return;
917 case GAL_WCS_COORDSYS_ECB1950:
918 *lng2p1=286.26975051; *lat2p1=38.34766498; *lng1p2=33.378919140; return;
919 case GAL_WCS_COORDSYS_ECJ2000:
920 *lng2p1=286.96654469; *lat2p1=38.34158720; *lng1p2=33.376119480; return;
921 case GAL_WCS_COORDSYS_GALACTIC:
922 *lng2p1=47.370000000; *lat2p1=6.320000000; *lng1p2=90.000000000; return;
923 case GAL_WCS_COORDSYS_SUPERGALACTIC:
924 *lng2p1=NAN; *lat2p1=NAN; *lng1p2=NAN; return;
925 default:
926 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
927 "fix the problem. The code '%d' isn't a recognized WCS "
928 "coordinate system ID for 'outsys' (input SUPERGALACTIC)", __func__,
929 PACKAGE_BUGREPORT, outsys);
930 }
931 break;
932 default:
933 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
934 "fix the problem. The code '%d' isn't a recognized WCS "
935 "coordinate system ID for 'insys'", __func__,
936 PACKAGE_BUGREPORT, insys);
937 }
938
939 }
940
941
942
943
944
945 static void
wcs_coordsys_ctypes(int coordsys,char ** clng,char ** clat,char ** radesys)946 wcs_coordsys_ctypes(int coordsys, char **clng, char **clat, char **radesys)
947 {
948 switch( coordsys)
949 {
950 case GAL_WCS_COORDSYS_EQB1950:
951 *clng="RA"; *clat="DEC"; *radesys="FK4"; break;
952 case GAL_WCS_COORDSYS_EQJ2000:
953 *clng="RA"; *clat="DEC"; *radesys="FK5"; break;
954 case GAL_WCS_COORDSYS_ECB1950:
955 *clng="ELON"; *clat="ELAT"; *radesys="FK4"; break;
956 case GAL_WCS_COORDSYS_ECJ2000:
957 *clng="ELON"; *clat="ELAT"; *radesys="FK5"; break;
958 case GAL_WCS_COORDSYS_GALACTIC:
959 *clng="GLON"; *clat="GLAT"; *radesys=NULL; break;
960 case GAL_WCS_COORDSYS_SUPERGALACTIC:
961 *clng="SLON"; *clat="SLAT"; *radesys=NULL; break;
962 default:
963 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
964 "fix the problem. The code '%d' isn't a recognized WCS "
965 "coordinate system ID for 'coordsys'", __func__,
966 PACKAGE_BUGREPORT, coordsys);
967 }
968 }
969
970
971
972
973 /* Convert the coordinate system. */
974 struct wcsprm *
gal_wcs_coordsys_convert(struct wcsprm * wcs,int outcoordsys)975 gal_wcs_coordsys_convert(struct wcsprm *wcs, int outcoordsys)
976 {
977 int incoordsys;
978 char *alt=NULL; /* Only concerned with primary wcs. */
979 double equinox=0.0f; /* To preserve current value. */
980 struct wcsprm *out=NULL;
981 double lng2p1=NAN, lat2p1=NAN, lng1p2=NAN;
982 char *clng=NULL, *clat=NULL, *radesys=NULL;
983
984 /* Just incase the input is a NULL pointer. */
985 if(wcs==NULL) return NULL;
986
987 /* Get the input's coordinate system and see if it should be converted at
988 all or not (if the output coordinate system is different). If its the
989 same, just copy the input and return. */
990 incoordsys=gal_wcs_coordsys_identify(wcs);
991 if(incoordsys==outcoordsys)
992 {
993 out=gal_wcs_copy(wcs);
994 return out;
995 }
996
997 /* Find the necessary pole coordinates. Note that we have already
998 accounted for the fact that the input and output coordinate systems
999 may be the same above, so the NaN outputs will never occur here. */
1000 wcs_coordsys_insys_pole_in_outsys(incoordsys, outcoordsys,
1001 &lng2p1, &lat2p1, &lng1p2);
1002
1003 /* Find the necessary CTYPE names of the output. */
1004 wcs_coordsys_ctypes(outcoordsys, &clng, &clat, &radesys);
1005
1006 /* Convert the WCS's coordinate system (if 'wcsccs' is available). */
1007 #if GAL_CONFIG_HAVE_WCSLIB_WCSCCS
1008 out=gal_wcs_copy(wcs);
1009 wcsccs(out, lng2p1, lat2p1, lng1p2, clng, clat, radesys, equinox, alt);
1010 #else
1011
1012 /* Just to avoid compiler warnings for 'equinox' and 'alt'. */
1013 if(alt) lng2p1+=equinox;
1014
1015 /* Print error message and abort. */
1016 error(EXIT_FAILURE, 0, "%s: the 'wcsccs' function isn't available "
1017 "in the version of WCSLIB that this Gnuastro was built with "
1018 "('wcsccs' was first available in WCSLIB 7.5, released on "
1019 "March 2021). Therefore, Gnuastro can't preform the WCS "
1020 "coordiante system conversion in the WCS. Please update your "
1021 "WCSLIB and re-build Gnuastro with it to use this feature. "
1022 "You can follow the instructions here to install the latest "
1023 "version of WCSLIB:\n"
1024 " https://www.gnu.org/software/gnuastro/manual/html_node/"
1025 "WCSLIB.html\n"
1026 "And then re-build Gnuastro as described here:\n"
1027 " https://www.gnu.org/software/gnuastro/manual/"
1028 "html_node/Quick-start.html\n\n",
1029 __func__);
1030 #endif
1031
1032 /* Return. */
1033 return out;
1034 }
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052 /*************************************************************
1053 *********** Distortions ***********
1054 *************************************************************/
1055 int
gal_wcs_distortion_from_string(char * distortion)1056 gal_wcs_distortion_from_string(char *distortion)
1057 {
1058 if( !strcmp(distortion,"TPD") ) return GAL_WCS_DISTORTION_TPD;
1059 else if( !strcmp(distortion,"SIP") ) return GAL_WCS_DISTORTION_SIP;
1060 else if( !strcmp(distortion,"TPV") ) return GAL_WCS_DISTORTION_TPV;
1061 else if( !strcmp(distortion,"DSS") ) return GAL_WCS_DISTORTION_DSS;
1062 else if( !strcmp(distortion,"WAT") ) return GAL_WCS_DISTORTION_WAT;
1063 else
1064 error(EXIT_FAILURE, 0, "WCS distortion name '%s' not recognized, "
1065 "currently recognized names are 'TPD', 'SIP', 'TPV', 'DSS' and "
1066 "'WAT'", distortion);
1067
1068 /* Control should not reach here. */
1069 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
1070 "problem. Control should not reach the end of this function",
1071 __func__, PACKAGE_BUGREPORT);
1072 return GAL_WCS_DISTORTION_INVALID;
1073 }
1074
1075
1076
1077
1078
1079 char *
gal_wcs_distortion_to_string(int distortion)1080 gal_wcs_distortion_to_string(int distortion)
1081 {
1082 /* Return the proper literal string. */
1083 switch(distortion)
1084 {
1085 case GAL_WCS_DISTORTION_TPD: return "TPD";
1086 case GAL_WCS_DISTORTION_SIP: return "SIP";
1087 case GAL_WCS_DISTORTION_TPV: return "TPV";
1088 case GAL_WCS_DISTORTION_DSS: return "DSS";
1089 case GAL_WCS_DISTORTION_WAT: return "WAT";
1090 default:
1091 error(EXIT_FAILURE, 0, "WCS distortion id '%d' isn't recognized",
1092 distortion);
1093 }
1094
1095 /* Control should not reach here. */
1096 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
1097 "problem. Control should not reach the end of this function",
1098 __func__, PACKAGE_BUGREPORT);
1099 return NULL;
1100 }
1101
1102
1103
1104
1105
1106 /* Check the type of distortion present and return the appropriate
1107 integer based on `enum gal_wcs_distortion`.
1108
1109 Parameters:
1110 struct wcsprm *wcs - The wcs parameters of the fits file.
1111
1112 Return:
1113 int out_distortion - The type of distortion present. */
1114 int
gal_wcs_distortion_identify(struct wcsprm * wcs)1115 gal_wcs_distortion_identify(struct wcsprm *wcs)
1116 {
1117 #if GAL_CONFIG_HAVE_WCSLIB_DIS_H
1118 struct disprm *dispre=NULL;
1119 struct disprm *disseq=NULL;
1120
1121 /* Small sanity check. */
1122 if(wcs==NULL) return GAL_WCS_DISTORTION_INVALID;
1123
1124 /* To help in reading. */
1125 disseq=wcs->lin.disseq;
1126 dispre=wcs->lin.dispre;
1127
1128 /* Check if distortion present. */
1129 if( disseq==NULL && dispre==NULL ) return GAL_WCS_DISTORTION_INVALID;
1130
1131 /* Check the type of distortion.
1132
1133 As mentioned in the WCS paper IV section 2.4.2 available at
1134 https://www.atnf.csiro.au/people/mcalabre/WCS/dcs_20040422.pdf, the
1135 DPja and DQia keywords are used to record the parameters required by
1136 the prior and sequent distortion functions respectively.
1137
1138 Now, as mentioned in dis.h file reference section in WCSLIB manual
1139 given here
1140 https://www.atnf.csiro.au/people/mcalabre/WCS/wcslib/dis_8h.html, TPV,
1141 DSS, and WAT are sequent polynomial distortions, while SIP is prior
1142 polynomial distortion. TPD is a superset of all these distortions and
1143 hence can be used both as a prior and sequent distortion polynomial.
1144
1145 References and citations:
1146 "Representations of distortions in FITS world coordinate systems",
1147 Calabretta, M.R. et al. (WCS Paper IV, draft dated 2004/04/22)
1148 */
1149
1150 if( dispre != NULL )
1151 {
1152 if( !strcmp(*dispre->dtype, "SIP") ) return GAL_WCS_DISTORTION_SIP;
1153 else if( !strcmp(*dispre->dtype, "TPD") ) return GAL_WCS_DISTORTION_TPD;
1154 else
1155 error(EXIT_FAILURE, 0, "%s: distortion '%s' isn't recognized in "
1156 "the 'dispre' structure of the given 'wcsprm'", __func__,
1157 *dispre->dtype);
1158 }
1159 else if( disseq != NULL )
1160 {
1161 if( !strcmp(*disseq->dtype, "TPV") ) return GAL_WCS_DISTORTION_TPV;
1162 else if( !strcmp(*disseq->dtype, "TPD") ) return GAL_WCS_DISTORTION_TPD;
1163 else if( !strcmp(*disseq->dtype, "DSS") ) return GAL_WCS_DISTORTION_DSS;
1164 else if( !strcmp(*disseq->dtype, "WAT") ) return GAL_WCS_DISTORTION_WAT;
1165 else
1166 error(EXIT_FAILURE, 0, "%s: distortion '%s' isn't recognized in "
1167 "the 'disseq' structure of the given 'wcsprm'", __func__,
1168 *dispre->dtype);
1169 }
1170
1171 /* Control should not reach here. */
1172 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to fix "
1173 "the problem. Control should not reach the end of this function",
1174 __func__, PACKAGE_BUGREPORT);
1175 #else
1176 /* The 'wcslib/dis.h' isn't present. */
1177 error(EXIT_FAILURE, 0, "%s: the installed version of WCSLIB on this "
1178 "system doesn't have the 'dis.h' header! Thus Gnuastro can't do "
1179 "distortion-related operations on the world coordinate system "
1180 "(WCS). To use these features, please upgrade your version "
1181 "of WCSLIB and re-build Gnuastro", __func__);
1182 #endif
1183 return GAL_WCS_DISTORTION_INVALID;
1184 }
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195 /* Convert a given distrotion type to other.
1196
1197 Parameters:
1198 struct wcsprm *wcs - The wcs parameters of the fits file.
1199 int out_distortion - The desired output distortion.
1200 size_t* fitsize - The size of the array along each dimension.
1201
1202 Return:
1203 struct wcsprm *outwcs - The transformed wcs parameters in the
1204 required distortion type. */
1205 struct wcsprm *
gal_wcs_distortion_convert(struct wcsprm * inwcs,int outdisptype,size_t * fitsize)1206 gal_wcs_distortion_convert(struct wcsprm *inwcs, int outdisptype,
1207 size_t *fitsize)
1208 {
1209 #if GAL_CONFIG_HAVE_WCSLIB_DIS_H
1210 struct wcsprm *outwcs=NULL;
1211 int indisptype=gal_wcs_distortion_identify(inwcs);
1212
1213 /* Make sure we have a PC+CDELT structure in the input WCS. */
1214 gal_wcs_decompose_pc_cdelt(inwcs);
1215
1216 /* If the input and output types are the same, just copy the input,
1217 otherwise, do the conversion. */
1218 if(indisptype==outdisptype) outwcs=gal_wcs_copy(inwcs);
1219 else
1220 switch(indisptype)
1221 {
1222 /* If there is no distortion in the input, just return a
1223 newly-allocated copy. */
1224 case GAL_WCS_DISTORTION_INVALID: outwcs=gal_wcs_copy(inwcs); break;
1225
1226 /* Input's distortion is SIP. */
1227 case GAL_WCS_DISTORTION_SIP:
1228 switch(outdisptype)
1229 {
1230 case GAL_WCS_DISTORTION_TPV:
1231 outwcs=gal_wcsdistortion_sip_to_tpv(inwcs);
1232 break;
1233 default:
1234 error(EXIT_FAILURE, 0, "%s: conversion from %s to %s is not yet "
1235 "supported. Please contact us at '%s'", __func__,
1236 gal_wcs_distortion_to_string(indisptype),
1237 gal_wcs_distortion_to_string(outdisptype),
1238 PACKAGE_BUGREPORT);
1239 }
1240 break;
1241
1242 /* Input's distortion is TPV. */
1243 case GAL_WCS_DISTORTION_TPV:
1244 switch(outdisptype)
1245 {
1246 case GAL_WCS_DISTORTION_SIP:
1247 if(fitsize==NULL)
1248 error(EXIT_FAILURE, 0, "%s: the size array is necessary "
1249 "for this conversion", __func__);
1250 outwcs=gal_wcsdistortion_tpv_to_sip(inwcs, fitsize);
1251 break;
1252 default:
1253 error(EXIT_FAILURE, 0, "%s: conversion from %s to %s is not yet "
1254 "supported. Please contact us at '%s'", __func__,
1255 gal_wcs_distortion_to_string(indisptype),
1256 gal_wcs_distortion_to_string(outdisptype),
1257 PACKAGE_BUGREPORT);
1258 }
1259 break;
1260
1261 /* Input's distortion is not yet supported.. */
1262 case GAL_WCS_DISTORTION_TPD:
1263 case GAL_WCS_DISTORTION_DSS:
1264 case GAL_WCS_DISTORTION_WAT:
1265 error(EXIT_FAILURE, 0, "%s: input %s distortion is not yet "
1266 "supported. Please contact us at '%s'", __func__,
1267 gal_wcs_distortion_to_string(indisptype),
1268 PACKAGE_BUGREPORT);
1269
1270 /* A bug! This distortion is not yet recognized. */
1271 default:
1272 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
1273 "the problem. The identifier '%d' is not recognized as a "
1274 "distortion", __func__, PACKAGE_BUGREPORT, indisptype);
1275 }
1276
1277 /* Return the converted WCS. */
1278 return outwcs;
1279 #else
1280 /* The 'wcslib/dis.h' isn't present. */
1281 error(EXIT_FAILURE, 0, "%s: the installed version of WCSLIB on this "
1282 "system doesn't have the 'dis.h' header! Thus Gnuastro can't do "
1283 "distortion-related operations on the world coordinate system "
1284 "(WCS). To use these features, please upgrade your version "
1285 "of WCSLIB and re-build Gnuastro", __func__);
1286 return NULL;
1287 #endif
1288 }
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310 /**************************************************************/
1311 /********** Utilities ************/
1312 /**************************************************************/
1313 /* Copy a given WSC structure into another one. */
1314 struct wcsprm *
gal_wcs_copy(struct wcsprm * wcs)1315 gal_wcs_copy(struct wcsprm *wcs)
1316 {
1317 struct wcsprm *out;
1318
1319 /* If the input WCS is NULL, return a NULL WCS. */
1320 if(wcs)
1321 {
1322 /* Allocate the output WCS structure. */
1323 errno=0;
1324 out=malloc(sizeof *out);
1325 if(out==NULL)
1326 error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for 'out'",
1327 __func__, sizeof *out);
1328
1329 /* Initialize the allocated WCS structure. The WCSLIB manual says "On
1330 the first invokation, and only the first invokation, wcsprm::flag
1331 must be set to -1 to initialize memory management"*/
1332 out->flag=-1;
1333 wcsini(1, wcs->naxis, out);
1334
1335 /* Copy the input WCS to the output WSC structure. */
1336 wcscopy(1, wcs, out);
1337 }
1338 else
1339 out=NULL;
1340
1341 /* Return the final output. */
1342 return out;
1343 }
1344
1345
1346
1347
1348
1349 /* Remove the algorithm part of CTYPE (anything after, and including, a
1350 '-') if necessary. */
1351 static void
wcs_ctype_noalgorithm(char * str)1352 wcs_ctype_noalgorithm(char *str)
1353 {
1354 size_t i, len=strlen(str);
1355 for(i=0;i<len;++i) if(str[i]=='-') { str[i]='\0'; break; }
1356 }
1357
1358
1359
1360
1361 /* See if the CTYPE string ends with TAN. */
1362 static int
wcs_ctype_has_tan(char * str)1363 wcs_ctype_has_tan(char *str)
1364 {
1365 size_t len=strlen(str);
1366
1367 return !strcmp(&str[len-3], "TAN");
1368 }
1369
1370
1371
1372
1373
1374 /* Remove dimension. */
1375 #define WCS_REMOVE_DIM_CHECK 0
1376 void
gal_wcs_remove_dimension(struct wcsprm * wcs,size_t fitsdim)1377 gal_wcs_remove_dimension(struct wcsprm *wcs, size_t fitsdim)
1378 {
1379 size_t c, i, j, naxis;
1380
1381 /* If the WCS structure is NULL, just return. */
1382 if(wcs==NULL) return;
1383
1384 /* Sanity check. */
1385 naxis=wcs->naxis;
1386 if(fitsdim==0 || fitsdim>naxis)
1387 error(EXIT_FAILURE, 0, "%s: requested dimension (fitsdim=%zu) must be "
1388 "larger than zero and smaller than the number of dimensions in "
1389 "the given WCS structure (%zu)", __func__, fitsdim, naxis);
1390
1391 /**************************************************/
1392 #if WCS_REMOVE_DIM_CHECK
1393 printf("\n\nfitsdim: %zu\n", fitsdim);
1394 printf("\n##################\n");
1395 /*
1396 wcs->pc[0]=0; wcs->pc[1]=1; wcs->pc[2]=2;
1397 wcs->pc[3]=3; wcs->pc[4]=4; wcs->pc[5]=5;
1398 wcs->pc[6]=6; wcs->pc[7]=7; wcs->pc[8]=8;
1399 */
1400 for(i=0;i<wcs->naxis;++i)
1401 {
1402 for(j=0;j<wcs->naxis;++j)
1403 printf("%-5g", wcs->pc[i*wcs->naxis+j]);
1404 printf("\n");
1405 }
1406 #endif
1407 /**************************************************/
1408
1409 /* First loop over the arrays. */
1410 for(i=0;i<naxis;++i)
1411 {
1412 /* The dimensions are in FITS order, but counting starts from 0, so
1413 we'll have to subtract 1 from 'fitsdim'. */
1414 if(i>fitsdim-1)
1415 {
1416 /* 1-D arrays. */
1417 if(wcs->crpix) wcs->crpix[i-1] = wcs->crpix[i];
1418 if(wcs->cdelt) wcs->cdelt[i-1] = wcs->cdelt[i];
1419 if(wcs->crval) wcs->crval[i-1] = wcs->crval[i];
1420 if(wcs->crota) wcs->crota[i-1] = wcs->crota[i];
1421 if(wcs->crder) wcs->crder[i-1] = wcs->crder[i];
1422 if(wcs->csyer) wcs->csyer[i-1] = wcs->csyer[i];
1423
1424 /* The strings are all statically allocated, so we don't need to
1425 check. */
1426 memcpy(wcs->cunit[i-1], wcs->cunit[i], 72);
1427 memcpy(wcs->ctype[i-1], wcs->ctype[i], 72);
1428 memcpy(wcs->cname[i-1], wcs->cname[i], 72);
1429
1430 /* For 2-D arrays, just bring up all the rows. We'll fix the
1431 columns in a second loop. */
1432 for(j=0;j<naxis;++j)
1433 {
1434 if(wcs->pc) wcs->pc[ (i-1)*naxis+j ] = wcs->pc[ i*naxis+j ];
1435 if(wcs->cd) wcs->cd[ (i-1)*naxis+j ] = wcs->cd[ i*naxis+j ];
1436 }
1437 }
1438 }
1439
1440
1441 /**************************************************/
1442 #if WCS_REMOVE_DIM_CHECK
1443 printf("\n###### Respective row removed (replaced).\n");
1444 for(i=0;i<wcs->naxis;++i)
1445 {
1446 for(j=0;j<wcs->naxis;++j)
1447 printf("%-5g", wcs->pc[i*wcs->naxis+j]);
1448 printf("\n");
1449 }
1450 #endif
1451 /**************************************************/
1452
1453
1454 /* Second loop for 2D arrays. */
1455 c=0;
1456 for(i=0;i<naxis;++i)
1457 for(j=0;j<naxis;++j)
1458 if(j!=fitsdim-1)
1459 {
1460 if(wcs->pc) wcs->pc[ c ] = wcs->pc[ i*naxis+j ];
1461 if(wcs->cd) wcs->cd[ c ] = wcs->cd[ i*naxis+j ];
1462 ++c;
1463 }
1464
1465
1466 /* Correct the total number of dimensions in the WCS structure. */
1467 naxis = wcs->naxis -= 1;
1468
1469
1470 /* The 'TAN' algorithm needs two dimensions. So we need to remove it when
1471 it can cause confusion. */
1472 switch(naxis)
1473 {
1474 /* The 'TAN' algorithm cannot be used for any single-dimensional
1475 dataset. So we'll have to remove it if it exists. */
1476 case 1:
1477 wcs_ctype_noalgorithm(wcs->ctype[0]);
1478 break;
1479
1480 /* For any other dimensionality, 'TAN' should be kept only when exactly
1481 two dimensions have it. */
1482 default:
1483
1484 c=0;
1485 for(i=0;i<naxis;++i)
1486 if( wcs_ctype_has_tan(wcs->ctype[i]) )
1487 ++c;
1488
1489 if(c!=2)
1490 for(i=0;i<naxis;++i)
1491 if( wcs_ctype_has_tan(wcs->ctype[i]) )
1492 wcs_ctype_noalgorithm(wcs->ctype[i]);
1493 break;
1494 }
1495
1496
1497
1498 /**************************************************/
1499 #if WCS_REMOVE_DIM_CHECK
1500 printf("\n###### Respective column removed.\n");
1501 for(i=0;i<naxis;++i)
1502 {
1503 for(j=0;j<naxis;++j)
1504 printf("%-5g", wcs->pc[i*naxis+j]);
1505 printf("\n");
1506 }
1507 printf("\n###### One final string\n");
1508 for(i=0;i<naxis;++i)
1509 printf("%s\n", wcs->ctype[i]);
1510 exit(0);
1511 #endif
1512 /**************************************************/
1513 }
1514
1515
1516
1517
1518
1519 /* Using the block data structure of the tile, add a WCS structure for
1520 it. In many cases, tiles are created for internal processing, so there
1521 is no need to keep their WCS. Hence for preformance reasons, when
1522 creating the tiles they don't have any WCS structure. When needed, this
1523 function can be used to add a WCS structure to the tile by copying the
1524 WCS structure of its block and correcting its starting points. If the
1525 tile already has a WCS structure, this function won't do anything.*/
1526 void
gal_wcs_on_tile(gal_data_t * tile)1527 gal_wcs_on_tile(gal_data_t *tile)
1528 {
1529 size_t i, start_ind, ndim=tile->ndim;
1530 gal_data_t *block=gal_tile_block(tile);
1531 size_t *coord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
1532 "coord");
1533
1534 /* If the tile already has a WCS structure, don't do anything. */
1535 if(tile->wcs) return;
1536 else
1537 {
1538 /* Copy the block's WCS into the tile. */
1539 tile->wcs=gal_wcs_copy(block->wcs);
1540
1541 /* Find the coordinates of the tile's starting index. */
1542 start_ind=gal_pointer_num_between(block->array, tile->array,
1543 block->type);
1544 gal_dimension_index_to_coord(start_ind, ndim, block->dsize, coord);
1545
1546 /* Correct the copied WCS structure. Note that crpix is indexed in
1547 the FITS/Fortran order while coord is ordered in C, it also starts
1548 counting from 1, not zero. */
1549 for(i=0;i<ndim;++i)
1550 tile->wcs->crpix[i] -= coord[ndim-1-i];
1551 /*
1552 printf("start_ind: %zu\n", start_ind);
1553 printf("coord: %zu, %zu\n", coord[1]+1, coord[0]+1);
1554 printf("CRPIX: %f, %f\n", tile->wcs->crpix[0], tile->wcs->crpix[1]);
1555 */
1556 }
1557
1558 /* Clean up. */
1559 free(coord);
1560 }
1561
1562
1563
1564
1565
1566 /* Return the Warping matrix of the given WCS structure. This will be the
1567 final matrix irrespective of the type of storage in the WCS
1568 structure. Recall that the FITS standard has several methods to store
1569 the matrix, which is up to this function to account for and return the
1570 final matrix. The output is an allocated DxD matrix where 'D' is the
1571 number of dimensions. */
1572 double *
gal_wcs_warp_matrix(struct wcsprm * wcs)1573 gal_wcs_warp_matrix(struct wcsprm *wcs)
1574 {
1575 double *out, crota2;
1576 size_t i, j, size=wcs->naxis*wcs->naxis;
1577
1578 /* Allocate the necessary array. */
1579 errno=0;
1580 out=malloc(size*sizeof *out);
1581 if(out==NULL)
1582 error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for 'out'",
1583 __func__, size*sizeof *out);
1584
1585 /* Fill in the array. */
1586 if(wcs->altlin & 0x1) /* Has a PCi_j array. */
1587 {
1588 for(i=0;i<wcs->naxis;++i)
1589 for(j=0;j<wcs->naxis;++j)
1590 out[i*wcs->naxis+j] = wcs->cdelt[i] * wcs->pc[i*wcs->naxis+j];
1591 }
1592 else if(wcs->altlin & 0x2) /* Has CDi_j array. */
1593 {
1594 for(i=0;i<size;++i)
1595 out[i]=wcs->cd[i];
1596 }
1597 else if(wcs->altlin & 0x4) /* Has CROTAi array. */
1598 {
1599 /* Basic sanity checks. */
1600 if(wcs->naxis!=2)
1601 error(EXIT_FAILURE, 0, "%s: CROTAi currently on works in 2 "
1602 "dimensions.", __func__);
1603 if(wcs->crota[0]!=0.0)
1604 error(EXIT_FAILURE, 0, "%s: CROTA1 is not zero", __func__);
1605
1606 /* CROTAi keywords are depreciated in the FITS standard. However, old
1607 files may still use them. For a full description of CROTAi
1608 keywords and their history (along with the conversion equations
1609 here), please see the link below:
1610
1611 https://fits.gsfc.nasa.gov/users_guide/users_guide/node57.html
1612
1613 Just note that the equations of the link above convert CROTAi to
1614 PC. But here we want the "final" matrix (after multiplication by
1615 the 'CDELT' values). So to speed things up, we won't bother
1616 dividing and then multiplying by the same CDELT values in the
1617 off-diagonal elements. */
1618 crota2=wcs->crota[1];
1619 out[0] = wcs->cdelt[0] * cos(crota2);
1620 out[1] = -1 * wcs->cdelt[1] *sin(crota2);
1621 out[2] = wcs->cdelt[0] * sin(crota2);
1622 out[3] = wcs->cdelt[1] * cos(crota2);
1623 }
1624 else
1625 error(EXIT_FAILURE, 0, "%s: currently only PCi_ja and CDi_ja keywords "
1626 "are recognized", __func__);
1627
1628 /* Return the result */
1629 return out;
1630 }
1631
1632
1633
1634
1635 /* Clean up small/negligible errros that are clearly caused by measurement
1636 errors in the PC and CDELT elements. */
1637 void
gal_wcs_clean_errors(struct wcsprm * wcs)1638 gal_wcs_clean_errors(struct wcsprm *wcs)
1639 {
1640 double crdcheck=NAN;
1641 size_t i, crdnum=0, ndim=wcs->naxis;
1642 double mean, crdsum=0, sum=0, min=FLT_MAX, max=0;
1643 double *pc=wcs->pc, *cdelt=wcs->cdelt, *crder=wcs->crder;
1644
1645 /* First clean up CDELT: if the CRDER keyword is set, then we'll use that
1646 as a reference, if not, we'll use the absolute floating point error
1647 defined in 'GAL_WCS_FLTERR'. */
1648 for(i=0; i<ndim; ++i)
1649 {
1650 sum+=cdelt[i];
1651 if(cdelt[i]>max) max=cdelt[i];
1652 if(cdelt[i]<min) min=cdelt[i];
1653 if(crder[i]!=UNDEFINED) {++crdnum; crdsum=crder[i];}
1654 }
1655 mean=sum/ndim;
1656 crdcheck = crdnum ? crdsum/crdnum : GAL_WCS_FLTERROR;
1657 if( (max-min)/mean < crdcheck )
1658 for(i=0; i<ndim; ++i)
1659 cdelt[i]=mean;
1660
1661 /* Now clean up the PC elements. If the diagonal elements are too close
1662 to 0, 1, or -1, set them to 0 or 1 or -1. */
1663 if(pc)
1664 for(i=0;i<ndim*ndim;++i)
1665 {
1666 if( fabs(pc[i] - 0 ) < GAL_WCS_FLTERROR ) pc[i]=0;
1667 else if( fabs(pc[i] - 1 ) < GAL_WCS_FLTERROR ) pc[i]=1;
1668 else if( fabs(pc[i] - -1 ) < GAL_WCS_FLTERROR ) pc[i]=-1;
1669 }
1670 }
1671
1672
1673
1674
1675
1676 /* According to the FITS standard, in the 'PCi_j' WCS formalism, the matrix
1677 elements m_{ij} are encoded in the 'PCi_j' keywords and the scale
1678 factors are encoded in the 'CDELTi' keywords. There is also another
1679 formalism (the 'CDi_j' formalism) which merges the two into one
1680 matrix.
1681
1682 However, WCSLIB's internal operations are apparently done in the 'PCi_j'
1683 formalism. So its outputs are also all in that format by default. When
1684 the input is a 'CDi_j', WCSLIB will still read the image into the
1685 'PCi_j' formalism and the 'CDELTi's are set to 1. This function will
1686 decompose the two matrices to give a reasonable 'CDELTi' and 'PCi_j' in
1687 such cases. */
1688 void
gal_wcs_decompose_pc_cdelt(struct wcsprm * wcs)1689 gal_wcs_decompose_pc_cdelt(struct wcsprm *wcs)
1690 {
1691 size_t i, j;
1692 double *ps, *warp;
1693
1694 /* If there is on WCS, then don't do anything. */
1695 if(wcs==NULL) return;
1696
1697 /* Get the pixel scale and full warp matrix. */
1698 warp=gal_wcs_warp_matrix(wcs);
1699 ps=gal_wcs_pixel_scale(wcs);
1700 if(ps==NULL) return;
1701
1702 /* For a check.
1703 printf("pc: %g, %g\n", pc[0], pc[1]);
1704 printf("warp: %g, %g, %g, %g\n", warp[0], warp[1],
1705 warp[2], warp[3]);
1706 */
1707
1708 /* Set the CDELTs. */
1709 for(i=0; i<wcs->naxis; ++i)
1710 wcs->cdelt[i] = ps[i];
1711
1712 /* Write the PC matrix. */
1713 for(i=0;i<wcs->naxis;++i)
1714 for(j=0;j<wcs->naxis;++j)
1715 wcs->pc[i*wcs->naxis+j] = warp[i*wcs->naxis+j]/ps[i];
1716
1717 /* According to the 'wcslib/wcs.h' header: "In particular, wcsset()
1718 resets wcsprm::cdelt to unity if CDi_ja is present (and no
1719 PCi_ja).". So apparently, when the input is a 'CDi_j', it might expect
1720 the 'CDELTi' elements to be 1.0. But we have changed that here, so we
1721 will correct the 'altlin' element of the WCS structure to make sure
1722 that WCSLIB only looks into the 'PCi_j' and 'CDELTi' and makes no
1723 assumptioins about 'CDELTi'. */
1724 wcs->altlin=1;
1725
1726 /* Clean up. */
1727 free(ps);
1728 free(warp);
1729 }
1730
1731
1732
1733
1734
1735 /* Set the WCS structure to use the CD matrix. */
1736 void
gal_wcs_to_cd(struct wcsprm * wcs)1737 gal_wcs_to_cd(struct wcsprm *wcs)
1738 {
1739 size_t i, j, naxis;
1740
1741 /* If there is on WCS, then don't do anything. */
1742 if(wcs==NULL) return;
1743
1744 /* 'wcs->altlin' identifies which rotation element is being used (PCi_j,
1745 CDi_J or CROTAi). For PCi_j, the first bit will be 1 (==1), for CDi_j,
1746 the second bit is 1 (==2) and for CROTAi, the third bit is 1 (==4). */
1747 naxis=wcs->naxis;
1748 switch(wcs->altlin)
1749 {
1750 /* PCi_j: Convert it to CDi_j. */
1751 case 1:
1752
1753 /* Fill in the CD matrix and correct the PC and CDELT arrays. We have
1754 to do this because ultimately, WCSLIB will be writing the PC and
1755 CDELT keywords, even when 'altlin' is 2. So effectively we have to
1756 multiply the PC and CDELT matrices, then set cdelt=1 in all
1757 dimensions. This is actually how WCSLIB reads a FITS header with
1758 only a CD matrix. */
1759 for(i=0;i<naxis;++i)
1760 {
1761 for(j=0;j<naxis;++j)
1762 wcs->cd[i*naxis+j] = wcs->pc[i*naxis+j] *= wcs->cdelt[i];
1763 wcs->cdelt[i]=1;
1764 }
1765
1766 /* Set the altlin to be the CD matrix and free the PC matrix. */
1767 wcs->altlin=2;
1768 break;
1769
1770 /* CDi_j: No need to do any conversion. */
1771 case 2: return; break;
1772
1773 /* CROTAi: not yet supported. */
1774 case 4:
1775 error(0, 0, "%s: WARNING: Conversion of 'CROTAi' keywords to the CD "
1776 "matrix is not yet supported (for lack of time!), please "
1777 "contact us at %s to add this feature. But this may not cause a "
1778 "problem at all, so please check if the output's WCS is "
1779 "reasonable", __func__, PACKAGE_BUGREPORT);
1780 break;
1781
1782 /* The value isn't supported! */
1783 default:
1784 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
1785 "problem. The value %d for wcs->altlin isn't recognized",
1786 __func__, PACKAGE_BUGREPORT, wcs->altlin);
1787 }
1788 }
1789
1790
1791
1792
1793
1794 /* The distance (along a great circle) on a sphere between two points
1795 is calculated here. Since the pixel sides are usually very small,
1796 we won't be using the direct formula:
1797
1798 cos(distance)=sin(d1)*sin(d2)+cos(d1)*cos(d2)*cos(r1-r2)
1799
1800 We will be using the haversine formula which better considering
1801 floating point errors (from Wikipedia:)
1802
1803 sin^2(distance)/2=sin^2( (d1-d2)/2 )+cos(d1)*cos(d2)*sin^2( (r1-r2)/2 )
1804
1805 Inputs and outputs are all in degrees.
1806 */
1807 double
gal_wcs_angular_distance_deg(double r1,double d1,double r2,double d2)1808 gal_wcs_angular_distance_deg(double r1, double d1, double r2, double d2)
1809 {
1810 /* Convert degrees to radians. */
1811 double r1r=r1*M_PI/180, d1r=d1*M_PI/180;
1812 double r2r=r2*M_PI/180, d2r=d2*M_PI/180;
1813
1814 /* To make things easier to read: */
1815 double a=sin( (d1r-d2r)/2 );
1816 double b=sin( (r1r-r2r)/2 );
1817
1818 /* Return the result: */
1819 return 2*asin( sqrt( a*a + cos(d1r)*cos(d2r)*b*b) ) * 180/M_PI;
1820 }
1821
1822
1823
1824
1825
1826 /* Return the pixel scale of the dataset in units of the WCS. */
1827 double *
gal_wcs_pixel_scale(struct wcsprm * wcs)1828 gal_wcs_pixel_scale(struct wcsprm *wcs)
1829 {
1830 gsl_vector S;
1831 gsl_matrix A, V;
1832 int warning_printed;
1833 gal_data_t *pixscale;
1834 size_t i, j, n, maxj, *permutation;
1835 double jvmax, *a, *out, *v, maxrow, minrow;
1836
1837 /* Only continue if a WCS exists. */
1838 if(wcs==NULL) return NULL;
1839
1840
1841 /* Write the full WCS rotation matrix into an array, irrespective of what
1842 style it was stored in the wcsprm structure ('PCi_j' style or 'CDi_j'
1843 style). */
1844 a=gal_wcs_warp_matrix(wcs);
1845
1846
1847 /* Now that everything is good, we can allocate the necessary memory. */
1848 n=wcs->naxis;
1849 v=gal_pointer_allocate(GAL_TYPE_FLOAT64, n*n, 0, __func__, "v");
1850 permutation=gal_pointer_allocate(GAL_TYPE_SIZE_T, n, 0, __func__,
1851 "permutation");
1852 pixscale=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &n, NULL,
1853 0, -1, 1, NULL, NULL, NULL);
1854
1855
1856 /* To avoid confusing issues with floating point errors being written in
1857 the non-diagonal elements of the FITS header PC or CD matrices, we
1858 need to check if the minimum and maximum values in each row are not
1859 several orders of magnitude apart.
1860
1861 Note that in some cases (for example a spectrum), one axis might be in
1862 degrees (value around 1e-5) and the other in angestroms (value around
1863 1e-10). So we can't look at the minimum and maximum of the whole
1864 matrix. However, in such cases, people will probably not warp/rotate
1865 the image to mix the coordinates. So the important thing to check is
1866 the minimum and maximum (non-zero) values in each row. */
1867 warning_printed=0;
1868 for(i=0;i<n;++i)
1869 {
1870 /* Find the minimum and maximum values in each row. */
1871 minrow=FLT_MAX;
1872 maxrow=-FLT_MAX;
1873 for(j=0;j<n;++j)
1874 if(a[i*n+j]!=0.0) /* We aren't concerned with 0 valued elements. */
1875 {
1876 /* We have to use absolutes because in cases like RA, the
1877 diagonal values in different rows may have different signs*/
1878 if(fabs(a[i*n+j])<minrow) minrow=fabs(a[i*n+j]);
1879 if(fabs(a[i*n+j])>maxrow) maxrow=fabs(a[i*n+j]);
1880 }
1881
1882 /* Do the check, print warning and make correction. */
1883 if(maxrow!=minrow
1884 && maxrow/minrow>1e5 /* The difference between elements is large */
1885 && maxrow/minrow<GAL_WCS_FLTERROR
1886 && warning_printed==0)
1887 {
1888 fprintf(stderr, "\nWARNING: The input WCS matrix (possibly taken "
1889 "from the FITS header keywords starting with 'CD' or 'PC') "
1890 "contains values with very different scales (more than "
1891 "10^5 different). This is probably due to floating point "
1892 "errors. These values might bias the pixel scale (and "
1893 "subsequent) calculations.\n\n"
1894 "You can see the respective matrix with one of the "
1895 "following two commands (depending on how the FITS file "
1896 "was written). Recall that if the desired extension/HDU "
1897 "isn't the default, you can choose it with the '--hdu' "
1898 "(or '-h') option before the '|' sign in these commands."
1899 "\n\n"
1900 " $ astfits file.fits -p | grep 'PC._.'\n"
1901 " $ astfits file.fits -p | grep 'CD._.'\n\n"
1902 "You can delete the ones with obvious floating point "
1903 "error values using the following command (assuming you "
1904 "want to delete 'CD1_2' and 'CD2_1'). Afterwards, you can "
1905 "re-run your original command to remove this warning "
1906 "message and possibly correct errors that it might have "
1907 "caused.\n\n"
1908 " $ astfits file.fits --delete=CD1_2 --delete=CD2_1\n\n"
1909 );
1910 warning_printed=1;
1911 }
1912 }
1913
1914
1915 /* Fill in the necessary GSL vector and Matrix structures. */
1916 S.size=n; S.stride=1; S.data=pixscale->array;
1917 V.size1=n; V.size2=n; V.tda=n; V.data=v;
1918 A.size1=n; A.size2=n; A.tda=n; A.data=a;
1919
1920
1921 /* Run GSL's Singular Value Decomposition, using one-sided Jacobi
1922 orthogonalization which computes the singular (scale) values to a
1923 higher relative accuracy. */
1924 gsl_linalg_SV_decomp_jacobi(&A, &V, &S);
1925
1926
1927 /* The raw pixel scale array produced from the singular value
1928 decomposition above is ordered based on values, not the input. So when
1929 the pixel scales in all the dimensions aren't the same (the units of
1930 the dimensions differ), the order of the values in 'pixelscale' will
1931 not necessarily correspond to the input's dimensions.
1932
1933 To correct the order, we can use the 'V' matrix to find the original
1934 position of the pixel scale values and then use permutation to
1935 re-order it correspondingly. The column with the largest (absolute)
1936 value will be taken as the one to be used for each row. */
1937 for(i=0;i<n;++i)
1938 {
1939 /* Find the column with the maximum value. */
1940 maxj=-1;
1941 jvmax=-FLT_MAX;
1942 for(j=0;j<n;++j)
1943 if(fabs(v[i*n+j])>jvmax)
1944 {
1945 maxj=j;
1946 jvmax=fabs(v[i*n+j]);
1947 }
1948
1949 /* Use the column with the maximum value for this dimension. */
1950 permutation[i]=maxj;
1951 }
1952
1953
1954 /* Apply the permutation described above. */
1955 gal_permutation_apply(pixscale, permutation);
1956
1957
1958 /* Clean up and return. */
1959 free(a);
1960 free(v);
1961 free(permutation);
1962 out=pixscale->array;
1963 pixscale->array=NULL;
1964 gal_data_free(pixscale);
1965 return out;
1966 }
1967
1968
1969
1970
1971
1972 /* Report the arcsec^2 area of the pixels in the image based on the
1973 WCS information in that image. */
1974 double
gal_wcs_pixel_area_arcsec2(struct wcsprm * wcs)1975 gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
1976 {
1977 double out;
1978 double *pixscale;
1979
1980 /* Some basic sanity checks. */
1981 if(wcs==NULL) return NAN;
1982 if(wcs->naxis==1) return NAN;
1983
1984 /* Check if the units of the axis are degrees or not. Currently all FITS
1985 images I have worked with use 'deg' for degrees. If other alternatives
1986 exist, we can add corrections later. */
1987 if( strcmp("deg", wcs->cunit[0]) || strcmp("deg", wcs->cunit[1]) )
1988 return NAN;
1989
1990 /* Get the pixel scales along each axis in degrees, then multiply. */
1991 pixscale=gal_wcs_pixel_scale(wcs);
1992 if(pixscale==NULL) return NAN;
1993
1994 /* Clean up and return the result. */
1995 out = pixscale[0] * pixscale[1] * 3600.0f * 3600.0f;
1996 free(pixscale);
1997 return out;
1998 }
1999
2000
2001
2002
2003
2004 int
gal_wcs_coverage(char * filename,char * hdu,size_t * ondim,double ** ocenter,double ** owidth,double ** omin,double ** omax)2005 gal_wcs_coverage(char *filename, char *hdu, size_t *ondim,
2006 double **ocenter, double **owidth, double **omin,
2007 double **omax)
2008 {
2009 fitsfile *fptr;
2010 struct wcsprm *wcs;
2011 int nwcs=0, type, status=0;
2012 char *name=NULL, *unit=NULL;
2013 gal_data_t *tmp, *coords=NULL;
2014 size_t i, ndim, *dsize=NULL, numrows;
2015 double *x=NULL, *y=NULL, *z=NULL, *min, *max, *center, *width;
2016
2017 /* Read the desired WCS (note that the linear matrix is irrelevant here,
2018 we'll just select PC because its the default WCS mode. */
2019 wcs=gal_wcs_read(filename, hdu, GAL_WCS_LINEAR_MATRIX_PC,
2020 0, 0, &nwcs);
2021
2022 /* If a WCS doesn't exist, return NULL. */
2023 if(wcs==NULL) return 0;
2024
2025 /* Make sure the input HDU is an image. */
2026 if( gal_fits_hdu_format(filename, hdu) != IMAGE_HDU )
2027 error(EXIT_FAILURE, 0, "%s (hdu %s): is not an image HDU, the "
2028 "'--skycoverage' option only applies to image extensions",
2029 filename, hdu);
2030
2031 /* Get the array information of the image. */
2032 fptr=gal_fits_hdu_open(filename, hdu, READONLY);
2033 gal_fits_img_info(fptr, &type, ondim, &dsize, &name, &unit);
2034 fits_close_file(fptr, &status);
2035 ndim=*ondim;
2036
2037 /* Abort if we have more than 3 dimensions. */
2038 if(ndim==1 || ndim>3) return 0;
2039
2040 /* Allocate the output datasets. */
2041 center=*ocenter=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim, 0, __func__,
2042 "ocenter");
2043 width=*owidth=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim, 0, __func__,
2044 "owidth");
2045 min=*omin=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim, 0, __func__,
2046 "omin");
2047 max=*omax=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim, 0, __func__,
2048 "omax");
2049
2050 /* Now that we have the number of dimensions in the image, allocate the
2051 space needed for the coordinates. */
2052 numrows = (ndim==2) ? 5 : 9;
2053 for(i=0;i<ndim;++i)
2054 {
2055 tmp=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numrows, NULL, 0,
2056 -1, 1, NULL, NULL, NULL);
2057 tmp->next=coords;
2058 coords=tmp;
2059 }
2060
2061 /* Fill in the coordinate arrays, Note that 'dsize' is ordered in C
2062 dimensions, for the WCS conversion, we need to have the dimensions
2063 ordered in FITS/Fortran order. */
2064 switch(ndim)
2065 {
2066 case 2:
2067 x=coords->array; y=coords->next->array;
2068 x[0] = 1; y[0] = 1;
2069 x[1] = dsize[1]; y[1] = 1;
2070 x[2] = 1; y[2] = dsize[0];
2071 x[3] = dsize[1]; y[3] = dsize[0];
2072 x[4] = dsize[1]/2 + (dsize[1]%2 ? 1 : 0.5f);
2073 y[4] = dsize[0]/2 + (dsize[0]%2 ? 1 : 0.5f);
2074 break;
2075 case 3:
2076 x=coords->array; y=coords->next->array; z=coords->next->next->array;
2077 x[0] = 1; y[0] = 1; z[0]=1;
2078 x[1] = dsize[2]; y[1] = 1; z[1]=1;
2079 x[2] = 1; y[2] = dsize[1]; z[2]=1;
2080 x[3] = dsize[2]; y[3] = dsize[1]; z[3]=1;
2081 x[4] = 1; y[4] = 1; z[4]=dsize[0];
2082 x[5] = dsize[2]; y[5] = 1; z[5]=dsize[0];
2083 x[6] = 1; y[6] = dsize[1]; z[6]=dsize[0];
2084 x[7] = dsize[2]; y[7] = dsize[1]; z[7]=dsize[0];
2085 x[8] = dsize[2]/2 + (dsize[2]%2 ? 1 : 0.5f);
2086 y[8] = dsize[1]/2 + (dsize[1]%2 ? 1 : 0.5f);
2087 z[8] = dsize[0]/2 + (dsize[0]%2 ? 1 : 0.5f);
2088 break;
2089 default:
2090 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
2091 "fix the problem. 'ndim' of %zu is not recognized.",
2092 __func__, PACKAGE_BUGREPORT, ndim);
2093 }
2094
2095 /* For a check:
2096 printf("IMAGE COORDINATES:\n");
2097 for(i=0;i<numrows;++i)
2098 if(ndim==2)
2099 printf("%-15g%-15g\n", x[i], y[i]);
2100 else
2101 printf("%-15g%-15g%-15g\n", x[i], y[i], z[i]);
2102 */
2103
2104 /* Convert to the world coordinate system. */
2105 gal_wcs_img_to_world(coords, wcs, 1);
2106
2107 /* For a check:
2108 printf("\nWORLD COORDINATES:\n");
2109 for(i=0;i<numrows;++i)
2110 if(ndim==2)
2111 printf("%-15g%-15g\n", x[i], y[i]);
2112 else
2113 printf("%-15g%-15g%-15g\n", x[i], y[i], z[i]);
2114 */
2115
2116 /* Get the minimum and maximum values in each dimension. */
2117 tmp=gal_statistics_minimum(coords);
2118 min[0] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
2119 tmp=gal_statistics_maximum(coords);
2120 max[0] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
2121 tmp=gal_statistics_minimum(coords->next);
2122 min[1] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
2123 tmp=gal_statistics_maximum(coords->next);
2124 max[1] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
2125 if(ndim>2)
2126 {
2127 tmp=gal_statistics_minimum(coords->next->next);
2128 min[2] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
2129 tmp=gal_statistics_maximum(coords->next->next);
2130 max[2] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
2131 }
2132
2133 /* Write the center and width. */
2134 switch(ndim)
2135 {
2136 case 2:
2137 center[0]=x[4]; center[1]=y[4];
2138 width[0]=max[0]-min[0]; width[1]=max[1]-min[1];
2139 break;
2140 case 3:
2141 center[0]=x[8]; center[1]=y[8]; center[2]=z[8];
2142 width[0]=max[0]-min[0]; width[1]=max[1]-min[1]; width[2]=max[2]-min[2];
2143 break;
2144 default:
2145 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to solve the "
2146 "problem. The value %zu is not a recognized dimension", __func__,
2147 PACKAGE_BUGREPORT, ndim);
2148 }
2149
2150 /* Clean up and return success. */
2151 free(dsize);
2152 wcsfree(wcs);
2153 return 1;
2154 }
2155
2156
2157
2158
2159
2160
2161
2162
2163
2164
2165
2166
2167
2168
2169
2170
2171
2172
2173
2174 /**************************************************************/
2175 /********** Array conversion ************/
2176 /**************************************************************/
2177 /* Some sanity checks for the WCS conversion functions. */
2178 static void
wcs_convert_sanity_check_alloc(gal_data_t * coords,struct wcsprm * wcs,const char * func,int ** stat,double ** phi,double ** theta,double ** world,double ** pixcrd,double ** imgcrd)2179 wcs_convert_sanity_check_alloc(gal_data_t *coords, struct wcsprm *wcs,
2180 const char *func, int **stat, double **phi,
2181 double **theta, double **world,
2182 double **pixcrd, double **imgcrd)
2183 {
2184 gal_data_t *tmp;
2185 size_t ndim=0, firstsize=0, size=coords->size;
2186
2187 /* Make sure a WCS structure is actually given. */
2188 if(wcs==NULL)
2189 error(EXIT_FAILURE, 0, "%s: input WCS structure is NULL", func);
2190
2191 for(tmp=coords; tmp!=NULL; tmp=tmp->next)
2192 {
2193 /* Count how many coordinates are given. */
2194 ++ndim;
2195
2196 /* Check the type of the input. */
2197 if(tmp->type!=GAL_TYPE_FLOAT64)
2198 error(EXIT_FAILURE, 0, "%s: input coordinates must have 'float64' "
2199 "type", func);
2200
2201 /* Make sure it has a single dimension. */
2202 if(tmp->ndim!=1)
2203 error(EXIT_FAILURE, 0, "%s: input coordinates for each dimension "
2204 "must each be one dimensional. Coordinate dataset %zu of the "
2205 "inputs has %zu dimensions", func, ndim, tmp->ndim);
2206
2207 /* See if all inputs have the same size. */
2208 if(ndim==1) firstsize=tmp->size;
2209 else
2210 if(firstsize!=tmp->size)
2211 error(EXIT_FAILURE, 0, "%s: all input coordinates must have the "
2212 "same number of elements. Coordinate dataset %zu has %zu "
2213 "elements while the first coordinate has %zu", func, ndim,
2214 tmp->size, firstsize);
2215 }
2216
2217 /* See if the number of coordinates given corresponds to the dimensions
2218 of the WCS structure. */
2219 if(ndim!=wcs->naxis)
2220 error(EXIT_FAILURE, 0, "%s: the number of input coordinates (%zu) does "
2221 "not match the dimensions of the input WCS structure (%d)", func,
2222 ndim, wcs->naxis);
2223
2224 /* Allocate all the necessary arrays. */
2225 *phi = gal_pointer_allocate( GAL_TYPE_FLOAT64, size, 0, __func__,
2226 "phi");
2227 *stat = gal_pointer_allocate( GAL_TYPE_INT32, size, 1, __func__,
2228 "stat");
2229 *theta = gal_pointer_allocate( GAL_TYPE_FLOAT64, size, 0, __func__,
2230 "theta");
2231 *world = gal_pointer_allocate( GAL_TYPE_FLOAT64, ndim*size, 0, __func__,
2232 "world");
2233 *imgcrd = gal_pointer_allocate( GAL_TYPE_FLOAT64, ndim*size, 0, __func__,
2234 "imgcrd");
2235 *pixcrd = gal_pointer_allocate( GAL_TYPE_FLOAT64, ndim*size, 0, __func__,
2236 "pixcrd");
2237 }
2238
2239
2240
2241
2242
2243 /* In Gnuastro, each column (coordinate for WCS conversion) is treated as a
2244 separate array in a 'gal_data_t' that are linked through a linked
2245 list. But in WCSLIB, the input is a single array (with multiple
2246 columns). This function will convert between the two. */
2247 static void
wcs_convert_list_to_from_array(gal_data_t * list,double * array,int * stat,size_t ndim,int to0from1)2248 wcs_convert_list_to_from_array(gal_data_t *list, double *array, int *stat,
2249 size_t ndim, int to0from1)
2250 {
2251 size_t i, d=0;
2252 gal_data_t *tmp;
2253
2254 for(tmp=list; tmp!=NULL; tmp=tmp->next)
2255 {
2256 /* Put all this coordinate's values into the single array that is
2257 input into or output from WCSLIB. */
2258 for(i=0;i<list->size;++i)
2259 {
2260 if(to0from1)
2261 ((double *)(tmp->array))[i] = stat[i] ? NAN : array[i*ndim+d];
2262 else
2263 array[i*ndim+d] = ((double *)(tmp->array))[i];
2264 }
2265
2266 /* Increment the dimension. */
2267 ++d;
2268 }
2269 }
2270
2271
2272
2273
2274
2275 /* Prepare the output of the WCS conversion functions. */
2276 static gal_data_t *
wcs_convert_prepare_out(gal_data_t * coords,struct wcsprm * wcs,int inplace)2277 wcs_convert_prepare_out(gal_data_t *coords, struct wcsprm *wcs, int inplace)
2278 {
2279 size_t i;
2280 gal_data_t *out=NULL;
2281 if(inplace)
2282 out=coords;
2283 else
2284 for(i=0;i<wcs->naxis;++i)
2285 gal_list_data_add_alloc(&out, NULL, GAL_TYPE_FLOAT64, 1,
2286 &coords->size, NULL, 0, coords->minmapsize,
2287 coords->quietmmap, wcs->ctype[i], wcs->cunit[i],
2288 NULL);
2289 return out;
2290 }
2291
2292
2293
2294
2295
2296 /* Convert world coordinates to image coordinates given the input WCS
2297 structure. The input must be a linked list of data structures of float64
2298 ('double') type. The top element of the linked list must be the first
2299 coordinate and etc. If 'inplace' is non-zero, then the output will be
2300 written into the input's allocated space. */
2301 gal_data_t *
gal_wcs_world_to_img(gal_data_t * coords,struct wcsprm * wcs,int inplace)2302 gal_wcs_world_to_img(gal_data_t *coords, struct wcsprm *wcs, int inplace)
2303 {
2304 gal_data_t *out;
2305 int status, *stat=NULL, ncoord=coords->size, nelem;
2306 double *phi=NULL, *theta=NULL, *world=NULL, *pixcrd=NULL, *imgcrd=NULL;
2307
2308 /* Some sanity checks. */
2309 wcs_convert_sanity_check_alloc(coords, wcs, __func__, &stat, &phi, &theta,
2310 &world, &pixcrd, &imgcrd);
2311 nelem=wcs->naxis; /* We have to make sure a WCS is given first. */
2312
2313
2314 /* Write the values from the input list of separate columns into a single
2315 array (WCSLIB input). */
2316 wcs_convert_list_to_from_array(coords, world, stat, wcs->naxis, 0);
2317
2318
2319 /* Use WCSLIB's wcss2p for the conversion. */
2320 status=wcss2p(wcs, ncoord, nelem, world, phi, theta, imgcrd, pixcrd, stat);
2321 if(status)
2322 error(EXIT_FAILURE, 0, "%s: wcss2p ERROR %d: %s", __func__, status,
2323 wcs_errmsg[status]);
2324
2325
2326 /* For a sanity check.
2327 {
2328 size_t i;
2329 printf("\n\n%s sanity check:\n", __func__);
2330 for(i=0;i<coords->size;++i)
2331 printf("(%g, %g, %g) --> (%g, %g, %g), [stat: %d]\n",
2332 world[i*3], world[i*3+1 ], world[i*3+2],
2333 pixcrd[i*3], pixcrd[i*3+1], pixcrd[i*3+2], stat[i]);
2334 }
2335 */
2336
2337
2338 /* Allocate the output arrays if they were not already allocated. */
2339 out=wcs_convert_prepare_out(coords, wcs, inplace);
2340
2341
2342 /* Write the output from a single array (WCSLIB output) into the output
2343 list of this function. */
2344 wcs_convert_list_to_from_array(out, pixcrd, stat, wcs->naxis, 1);
2345
2346
2347 /* Clean up. */
2348 free(phi);
2349 free(stat);
2350 free(theta);
2351 free(world);
2352 free(imgcrd);
2353 free(pixcrd);
2354
2355 /* Return the output list of coordinates. */
2356 return out;
2357 }
2358
2359
2360
2361
2362
2363 /* Similar to 'gal_wcs_world_to_img'. */
2364 gal_data_t *
gal_wcs_img_to_world(gal_data_t * coords,struct wcsprm * wcs,int inplace)2365 gal_wcs_img_to_world(gal_data_t *coords, struct wcsprm *wcs, int inplace)
2366 {
2367 gal_data_t *out;
2368 int status, *stat=NULL, ncoord=coords->size, nelem;
2369 double *phi=NULL, *theta=NULL, *world=NULL, *pixcrd=NULL, *imgcrd=NULL;
2370
2371 /* Some sanity checks. */
2372 wcs_convert_sanity_check_alloc(coords, wcs, __func__, &stat, &phi, &theta,
2373 &world, &pixcrd, &imgcrd);
2374 nelem=wcs->naxis; /* We have to make sure a WCS is given first. */
2375
2376
2377 /* Write the values from the input list of separate columns into a single
2378 array (WCSLIB input). */
2379 wcs_convert_list_to_from_array(coords, pixcrd, stat, wcs->naxis, 0);
2380
2381
2382 /* Use WCSLIB's wcsp2s for the conversion. */
2383 status=wcsp2s(wcs, ncoord, nelem, pixcrd, imgcrd, phi, theta, world, stat);
2384 if(status)
2385 error(EXIT_FAILURE, 0, "%s: wcsp2s ERROR %d: %s", __func__, status,
2386 wcs_errmsg[status]);
2387
2388
2389 /* For a check.
2390 {
2391 size_t i;
2392 printf("\n\n%s sanity check (%d dimensions):\n", __func__, nelem);
2393 for(i=0;i<coords->size;++i)
2394 switch(nelem)
2395 {
2396 case 2:
2397 printf("(%-10g %-10g) --> (%-10g %-10g), [stat: %d]\n",
2398 pixcrd[i*2], pixcrd[i*2+1],
2399 world[i*2], world[i*2+1], stat[i]);
2400 break;
2401 case 3:
2402 printf("(%g, %g, %g) --> (%g, %g, %g), [stat: %d]\n",
2403 pixcrd[i*3], pixcrd[i*3+1], pixcrd[i*3+2],
2404 world[i*3], world[i*3+1], world[i*3+2], stat[i]);
2405 break;
2406 }
2407 }
2408 */
2409
2410
2411 /* Allocate the output arrays if they were not already allocated. */
2412 out=wcs_convert_prepare_out(coords, wcs, inplace);
2413
2414
2415 /* Write the output from a single array (WCSLIB output) into the output
2416 list of this function. */
2417 wcs_convert_list_to_from_array(out, world, stat, wcs->naxis, 1);
2418
2419
2420 /* Clean up. */
2421 free(phi);
2422 free(stat);
2423 free(theta);
2424 free(world);
2425 free(imgcrd);
2426 free(pixcrd);
2427
2428
2429 /* Return the output list of coordinates. */
2430 return out;
2431 }
2432