1 /*********************************************************************
2 Warp - Warp images using projective mapping.
3 Warp 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) 2016-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 <argp.h>
26 #include <errno.h>
27 #include <error.h>
28 #include <stdio.h>
29 #include <string.h>
30
31 #include <gnuastro/wcs.h>
32 #include <gnuastro/fits.h>
33 #include <gnuastro/array.h>
34 #include <gnuastro/table.h>
35 #include <gnuastro/threads.h>
36 #include <gnuastro/pointer.h>
37
38 #include <gnuastro-internal/timing.h>
39 #include <gnuastro-internal/options.h>
40 #include <gnuastro-internal/checkset.h>
41 #include <gnuastro-internal/fixedstringmacros.h>
42
43 #include "main.h"
44
45 #include "ui.h"
46 #include "authors-cite.h"
47
48
49
50
51
52 /**************************************************************/
53 /********* Argp necessary global entities ************/
54 /**************************************************************/
55 /* Definition parameters for the Argp: */
56 const char *
57 argp_program_version = PROGRAM_STRING "\n"
58 GAL_STRINGS_COPYRIGHT
59 "\n\nWritten/developed by "PROGRAM_AUTHORS;
60
61 const char *
62 argp_program_bug_address = PACKAGE_BUGREPORT;
63
64 static char
65 args_doc[] = "ASTRdata";
66
67 const char
68 doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will warp/transform the "
69 "input image using an input coordinate matrix. Currently it accepts any "
70 "general projective mapping (which includes affine mappings as a "
71 "subset). \n"
72 GAL_STRINGS_MORE_HELP_INFO
73 /* After the list of options: */
74 "\v"
75 PACKAGE_NAME" home page: "PACKAGE_URL;
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96 /**************************************************************/
97 /********* Initialize & Parse command-line **************/
98 /**************************************************************/
99 static void
ui_initialize_options(struct warpparams * p,struct argp_option * program_options,struct argp_option * gal_commonopts_options)100 ui_initialize_options(struct warpparams *p,
101 struct argp_option *program_options,
102 struct argp_option *gal_commonopts_options)
103 {
104 size_t i;
105 struct gal_options_common_params *cp=&p->cp;
106
107
108 /* Set the necessary common parameters structure. */
109 cp->program_struct = p;
110 cp->program_name = PROGRAM_NAME;
111 cp->program_exec = PROGRAM_EXEC;
112 cp->program_bibtex = PROGRAM_BIBTEX;
113 cp->program_authors = PROGRAM_AUTHORS;
114 cp->poptions = program_options;
115 cp->numthreads = gal_threads_number();
116 cp->coptions = gal_commonopts_options;
117
118
119 /* Set the mandatory common options. */
120 for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
121 {
122 /* Select individually. */
123 switch(cp->coptions[i].key)
124 {
125 case GAL_OPTIONS_KEY_MINMAPSIZE:
126 cp->coptions[i].mandatory=GAL_OPTIONS_MANDATORY;
127 break;
128
129 case GAL_OPTIONS_KEY_SEARCHIN:
130 case GAL_OPTIONS_KEY_TABLEFORMAT:
131 case GAL_OPTIONS_KEY_STDINTIMEOUT:
132 cp->coptions[i].flags=OPTION_HIDDEN;
133 break;
134 }
135
136 /* Select by group. */
137 switch(cp->coptions[i].group)
138 {
139 case GAL_OPTIONS_GROUP_TESSELLATION:
140 cp->coptions[i].doc=NULL; /* Necessary to remove title. */
141 cp->coptions[i].flags=OPTION_HIDDEN;
142 break;
143 }
144 }
145 }
146
147
148
149
150
151 /* Parse a single option: */
152 error_t
parse_opt(int key,char * arg,struct argp_state * state)153 parse_opt(int key, char *arg, struct argp_state *state)
154 {
155 struct warpparams *p = state->input;
156
157 /* Pass 'gal_options_common_params' into the child parser. */
158 state->child_inputs[0] = &p->cp;
159
160 /* In case the user incorrectly uses the equal sign (for example
161 with a short format or with space in the long format, then 'arg'
162 start with (if the short version was called) or be (if the long
163 version was called with a space) the equal sign. So, here we
164 check if the first character of arg is the equal sign, then the
165 user is warned and the program is stopped: */
166 if(arg && arg[0]=='=')
167 argp_error(state, "incorrect use of the equal sign ('='). For short "
168 "options, '=' should not be used and for long options, "
169 "there should be no space between the option, equal sign "
170 "and value");
171
172 /* Set the key to this option. */
173 switch(key)
174 {
175 /* Read the non-option tokens (arguments): */
176 case ARGP_KEY_ARG:
177 /* The user may give a shell variable that is empty! In that case
178 'arg' will be an empty string! We don't want to account for such
179 cases (and give a clear error that no input has been given). */
180 if(p->inputname)
181 argp_error(state, "only one argument (input file) should be given");
182 else
183 if(arg[0]!='\0') p->inputname=arg;
184 break;
185
186 /* This is an option, set its value. */
187 default:
188 return gal_options_set_from_key(key, arg, p->cp.poptions, &p->cp);
189 }
190
191 return 0;
192 }
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212 /**************************************************************/
213 /********** Modular matrix linked list *************/
214 /**************************************************************/
215 /* Save the codes of the user's desired modular warpings into the linked
216 list. Because the types of these options are 'GAL_TYPE_INVALID', this
217 function will not be called when printing the full list of parameters
218 and their values. */
219 static void *
ui_add_to_modular_warps_ll(struct argp_option * option,char * arg,char * filename,size_t lineno,void * params)220 ui_add_to_modular_warps_ll(struct argp_option *option, char *arg,
221 char *filename, size_t lineno, void *params)
222 {
223 size_t i;
224 double tmp;
225 gal_data_t *new;
226 struct warpparams *p=(struct warpparams *)params;
227
228 /* When an argument is necessary (note that '--align' doesn't take an
229 argument), make sure we actually have a string to parse. Also note
230 that if an argument is necessary, but none is given Argp will
231 automatically abort the program with an informative error. */
232 if(arg && *arg=='\0')
233 error(EXIT_FAILURE, 0, "empty string given to '--%s'", option->name);
234
235 /* Parse the (possible) arguments. */
236 if(option->key==UI_KEY_ALIGN)
237 {
238 /* For functions the standard checking isn't done, so first, we'll
239 make sure that if we are in a configuration file (where
240 'arg!=NULL'), the value is either 0 or 1. */
241 if( arg && strcmp(arg, "0") && strcmp(arg, "1") )
242 error_at_line(EXIT_FAILURE, 0, filename, lineno, "the '--align' "
243 "option takes no arguments. In a configuration file "
244 "it can only have the values '1' or '0', indicating "
245 "if it should be used or not");
246
247 /* Align doesn't take any values, but if called in a configuration
248 file with a value of '0', we should ignore it. */
249 if(arg && *arg=='0') return NULL;
250
251 /* Allocate the data structure. */
252 new=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 0, NULL, NULL, 0, -1, 1,
253 NULL, NULL, NULL);
254 }
255 else new=gal_options_parse_list_of_numbers(arg, filename, lineno);
256
257
258 /* If this was a matrix, then put it in the matrix element of the main
259 data structure. Otherwise, add the list of given values to the modular
260 warpings list. */
261 if(option->key==UI_KEY_MATRIX)
262 {
263 /* Some sanity checks. */
264 if(p->matrix)
265 error_at_line(EXIT_FAILURE, 0, filename, lineno, "only one matrix "
266 "may be given, you can use multiple modular warpings");
267 if(new->size!=4 && new->size!=9)
268 error_at_line(EXIT_FAILURE, 0, filename, lineno, "only a 4 or 9 "
269 "element 'matrix' is currently acceptable. '%s' has "
270 "%zu elements", arg, new->size);
271
272 /* Keep the matrix in the main structure. */
273 p->matrix=new;
274 }
275 else
276 {
277 /* No more than two numbers should be given for the modular
278 warpings. */
279 if(new->size>2)
280 error_at_line(EXIT_FAILURE, 0, filename, lineno, "%zu numbers "
281 "given to the '%s' option. Modular warpings can "
282 "accept 2 numbers at the most currently (for 2D "
283 "datasets)", new->size, option->name);
284
285 /* Some modular-warp specific sanity checks: rotate only needs one
286 number, and flip's values should only be 0 and 1. */
287 if(option->key==UI_KEY_ROTATE)
288 {
289 if(new->size!=1)
290 error_at_line(EXIT_FAILURE, 0, filename, lineno, "the 'rotate' "
291 "option only takes one value (the angle of rotation). "
292 "You have given: '%s'", arg);
293 }
294 else if (option->key==UI_KEY_FLIP)
295 {
296 for(i=0;i<new->size;++i)
297 {
298 tmp=((double *)(new->array))[i];
299 if(tmp!=0.0f && tmp!=1.0f)
300 error_at_line(EXIT_FAILURE, 0, filename, lineno, "'flip' "
301 "only takes values of '1' and '0'. You have "
302 "given '%s'", arg);
303 }
304 }
305
306 /* Keep the final value. */
307 new->status=option->key;
308 new->next=p->modularll;
309 p->modularll=new;
310 }
311 return NULL;
312 }
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333 /**************************************************************/
334 /*************** Sanity Check *******************/
335 /**************************************************************/
336 static void
ui_check_options_and_arguments(struct warpparams * p)337 ui_check_options_and_arguments(struct warpparams *p)
338 {
339 /* Read the input.*/
340 if(p->inputname)
341 {
342 /* Make sure a HDU is given. */
343 if( gal_fits_file_recognized(p->inputname) && p->cp.hdu==NULL )
344 error(EXIT_FAILURE, 0, "no HDU specified, you can use the '--hdu' "
345 "('-h') option and give it the HDU number (starting from "
346 "zero), or extension name (generally, anything acceptable "
347 "by CFITSIO)");
348
349 /* Read the input image as double type and its WCS structure. */
350 p->input=gal_array_read_one_ch_to_type(p->inputname, p->cp.hdu,
351 NULL, GAL_TYPE_FLOAT64,
352 p->cp.minmapsize,
353 p->cp.quietmmap);
354
355 /* Read the WCS and remove one-element wide dimension(s). */
356 p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu,
357 p->cp.wcslinearmatrix, p->hstartwcs,
358 p->hendwcs, &p->input->nwcs);
359 p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
360 p->input->dsize,
361 p->input->wcs);
362
363 /* Currently Warp only works on 2D images. */
364 if(p->input->ndim!=2)
365 error(EXIT_FAILURE, 0, "input has %zu dimensions but Warp currently "
366 "only works on 2D datasets (images).\n\n"
367 "We do plan to add 3D functionality (see "
368 "https://savannah.gnu.org/task/?15729), so please get in "
369 "touch if you need it (any further interest, support or help "
370 "would be useful)", p->input->ndim);
371
372 /* Get basic WCS information. */
373 if(p->input->wcs)
374 {
375 p->pixelscale=gal_wcs_pixel_scale(p->input->wcs);
376 if(p->pixelscale==NULL)
377 error(EXIT_FAILURE, 0, "%s (hdu %s): the pixel scale couldn't "
378 "be deduced from the WCS.", p->inputname, p->cp.hdu);
379 p->inwcsmatrix=gal_wcs_warp_matrix(p->input->wcs);
380 }
381 }
382 else
383 error(EXIT_FAILURE, 0, "no input file is specified");
384 }
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405 /**************************************************************/
406 /*************** Matrix preparations ******************/
407 /**************************************************************/
408 static void
ui_error_no_warps()409 ui_error_no_warps()
410 {
411 error(EXIT_FAILURE, 0, "no warping specified, you can either use the "
412 "'--matrix' option for any low-level warp, or specify multipole "
413 "modular warpings with options like '--rotate', '--scale' and etc. "
414 "You can see the full list with the '--help' option");
415 }
416
417
418
419
420
421 /* This function is mainly for easy checking/debugging. */
422 static void
ui_matrix_print(double * matrix)423 ui_matrix_print(double *matrix)
424 {
425 printf("%-10.3f%-10.3f%-10.3f\n", matrix[0], matrix[1], matrix[2]);
426 printf("%-10.3f%-10.3f%-10.3f\n", matrix[3], matrix[4], matrix[5]);
427 printf("%-10.3f%-10.3f%-10.3f\n", matrix[6], matrix[7], matrix[8]);
428 }
429
430
431
432
433
434 static void
ui_matrix_prepare_raw(struct warpparams * p)435 ui_matrix_prepare_raw(struct warpparams *p)
436 {
437 size_t *dsize;
438 double *in=p->matrix->array, *final;
439
440 /* If the matrix was 2D, then convert it to 3D. Note that we done a size
441 check when reading the matrix, so at this point, it either has 9
442 elements, or 4. */
443 if(p->matrix->size==4)
444 {
445 /* Allocate the final matrix. */
446 final=gal_pointer_allocate(GAL_TYPE_FLOAT64, 9, 0, __func__, "final");
447
448 /* Fill in the final 3x3 matrix from the 2x2 matrix. */
449 final[0]=in[0]; final[1]=in[1]; final[2]=0.0f;
450 final[3]=in[2]; final[4]=in[3]; final[5]=0.0f;
451 final[6]=0.0f; final[7]=0.0f; final[8]=1.0f;
452
453 /* Free the old matrix array and put in the new one. */
454 free(p->matrix->array);
455 p->matrix->size=9;
456 p->matrix->array=final;
457 }
458
459 /* Correct the dimensional information, because the matrix was read as a
460 single dimensional list of numbers. */
461 free(p->matrix->dsize);
462 dsize=p->matrix->dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, 2, 0,
463 __func__, "dsize");
464 dsize[0]=dsize[1]=3;
465 p->matrix->ndim=2;
466 }
467
468
469
470
471
472 /* Set the matrix so the image is aligned with the axises. Note that
473 WCSLIB automatically fills the CRPI */
474 static void
ui_matrix_make_align(struct warpparams * p,double * tmatrix)475 ui_matrix_make_align(struct warpparams *p, double *tmatrix)
476 {
477 double A, *w, *P, x[4];
478
479 /* Make sure the input image had a WCS structure. */
480 if(p->input->wcs==NULL)
481 error(EXIT_FAILURE, 0, "%s (hdu: %s): no WCS information present, "
482 "hence the '--align' option cannot be used", p->inputname,
483 p->cp.hdu);
484
485 /* Check if there is only two WCS axises: */
486 if(p->input->wcs->naxis!=2)
487 error(EXIT_FAILURE, 0, "the WCS structure of %s (hdu: %s) has %d "
488 "axises. For the '--align' option to operate it must be 2",
489 p->inputname, p->cp.hdu, p->input->wcs->naxis);
490
491
492 /* For help in reading, use aliases for the WCS matrix and pixel scale.*/
493 P=p->pixelscale;
494 w=p->inwcsmatrix;
495
496
497 /* Lets call the given WCS orientation 'W', the rotation matrix we want
498 to find as 'X' and the final (aligned matrix) 'P' (which is the pixel
499 scale):
500
501 X W P
502 ------ ------ -------
503 x0 x1 w0 w1 -P0 0
504 x2 x3 * w2 w3 = 0 P1
505
506 Let's open up the matrix multiplication, so we can find the 'X'
507 elements as function of the 'W' elements and 'a'.
508
509 x0*w0 + x1*w2 = -P0 (1)
510 x0*w1 + x1*w3 = 0 (2)
511 x2*w0 + x3*w2 = 0 (3)
512 x2*w1 + x3*w3 = P1 (4)
513
514 Let's bring the X with the smaller index in each equation to the left
515 side:
516
517 x0 = (-w2/w0)*x1 - P0/w0 (5)
518 x0 = (-w3/w1)*x1 (6)
519 x2 = (-w2/w0)*x3 (7)
520 x2 = (-w3/w1)*x3 + P1/w1 (8)
521
522 Using (5) and (6) we can find x0 and x1, by first eliminating x0:
523
524 (-w2/w0)*x1 - P0/w0 = (-w3/w1)*x1 -> (w3/w1 - w2/w0) * x1 = P0/w0
525
526 For easy reading/writing, let's define: A = (w3/w1 - w2/w0)
527
528 --> x1 = P0 / w0 / A
529 --> x0 = -1 * x1 * w3 / w1
530
531 Similar to the above, we can find x2 and x3 from (7) and (8):
532
533 (-w2/w0)*x3 = (-w3/w1)*x3 + P1/w1 -> (w3/w1 - w2/w0) * x3 = P1/w1
534
535 --> x3 = P1 / w1 / A
536 --> x2 = -1 * x3 * w2 / w0
537
538 Note that when the image is already aligned (off-diagonals are zero),
539 only the signs of the diagonal elements matter. */
540 if( w[1]==0.0f && w[2]==0.0f )
541 {
542 x[0] = w[0]<0 ? 1.0f : -1.0f; /* Has to be negative. */
543 x[1] = 0.0f;
544 x[2] = 0.0f;
545 x[3] = w[3]>0 ? 1.0f : -1.0f; /* Has to be positive. */
546 }
547 else if (w[0]==0.0f && w[3]==0.0f )
548 {
549 x[0] = 0.0f;
550 x[1] = w[1]<0 ? 1.0f : -1.0f; /* Has to be negative. */
551 x[2] = w[2]>0 ? 1.0f : -1.0f; /* Has to be positive. */
552 x[3] = 0.0f;
553 }
554 else
555 {
556 A = (w[3]/w[1]) - (w[2]/w[0]);
557 x[1] = P[0] / w[0] / A;
558 x[3] = P[1] / w[1] / A;
559 x[0] = -1 * x[1] * w[3] / w[1];
560 x[2] = -1 * x[3] * w[2] / w[0];
561
562 /* When the input matrix is flipped, the diagonal elements of the
563 necessary rotation will have a different sign. */
564 if(x[0]*x[3]<0)
565 if(!p->cp.quiet)
566 error(0, 0, "%s (hdu %s): WARNING (bug #55217)! The alignment may "
567 "not be correct! This is a recognized bug, we just haven't "
568 "had enough time to fix it yet, any help or support is "
569 "most welcome.\n\n"
570 "This may be due to the East and North not being left-hand "
571 "oriented. If this is the case, you can fix it by flipping "
572 "the image along the problematic axis (with '--flip=1,0' or "
573 "'--flip=0,1') and re-running Warp with '--align'. Note "
574 "that you can't yet call '--flip' in the same command as "
575 "'--align'", p->inputname, p->cp.hdu);
576 }
577
578 /* For a check:
579 printf("ps: (%e, %e)\n", P[0], P[1]);
580 printf("w:\n");
581 printf(" %.8e %.8e\n", w[0], w[1]);
582 printf(" %.8e %.8e\n", w[2], w[3]);
583 printf("x:\n");
584 printf(" %.8e %.8e\n", x[0], x[1]);
585 printf(" %.8e %.8e\n", x[2], x[3]);
586 exit(0);
587 */
588
589 /* Put the matrix elements into the output array: */
590 tmatrix[0]=x[0]; tmatrix[1]=x[1]; tmatrix[2]=0.0f;
591 tmatrix[3]=x[2]; tmatrix[4]=x[3]; tmatrix[5]=0.0f;
592 tmatrix[6]=0.0f; tmatrix[7]=0.0f; tmatrix[8]=1.0f;
593 }
594
595
596
597
598
599 static void
ui_matrix_inplacw_multiply(double * in,double * with)600 ui_matrix_inplacw_multiply(double *in, double *with)
601 {
602 /* 'tin' will keep the values of the input array because we want to
603 write the multiplication result in the input array. */
604 double tin[9]={in[0],in[1],in[2],in[3],in[4],in[5],in[6],in[7],in[8]};
605
606 /* For easy checking, here are the matrix/memory layouts:
607 tin[0] tin[1] tin[2] with[0] with[1] with[2]
608 tin[3] tin[4] tin[5] * with[3] with[4] with[5]
609 tin[6] tin[7] tin[8] with[6] with[7] with[8] */
610 in[0] = tin[0]*with[0] + tin[1]*with[3] + tin[2]*with[6];
611 in[1] = tin[0]*with[1] + tin[1]*with[4] + tin[2]*with[7];
612 in[2] = tin[0]*with[2] + tin[1]*with[5] + tin[2]*with[8];
613
614 in[3] = tin[3]*with[0] + tin[4]*with[3] + tin[5]*with[6];
615 in[4] = tin[3]*with[1] + tin[4]*with[4] + tin[5]*with[7];
616 in[5] = tin[3]*with[2] + tin[4]*with[5] + tin[5]*with[8];
617
618 in[6] = tin[6]*with[0] + tin[7]*with[3] + tin[8]*with[6];
619 in[7] = tin[6]*with[1] + tin[7]*with[4] + tin[8]*with[7];
620 in[8] = tin[6]*with[2] + tin[7]*with[5] + tin[8]*with[8];
621 }
622
623
624
625
626
627
628 static void
ui_matrix_from_modular(struct warpparams * p)629 ui_matrix_from_modular(struct warpparams *p)
630 {
631 gal_data_t *pop;
632 size_t dsize[]={3,3};
633 double s, c, v1, v2, *final, module[9]={1,0,0, 0,1,0, 0,0,1};
634
635 /* Reverse the list of modular warpings to be in the same order as the
636 user specified.*/
637 gal_list_data_reverse(&p->modularll);
638
639 /* Allocate space for the final matrix. */
640 p->matrix=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 2, dsize, NULL, 0,
641 p->cp.minmapsize, p->cp.quietmmap,
642 NULL, NULL, NULL);
643 final=p->matrix->array;
644
645 /* Fill in the final matrix to start with. */
646 final[0]=1.0f; final[1]=0.0f; final[2]=0.0f;
647 final[3]=0.0f; final[4]=1.0f; final[5]=0.0f;
648 final[6]=0.0f; final[7]=0.0f; final[8]=1.0f;
649
650 /* Apply all modular warps. */
651 while(p->modularll)
652 {
653 /* Pop the top element. */
654 pop=gal_list_data_pop(&p->modularll);
655
656 /* Set the (possibly) two values given for this warp. */
657 v1 = pop->ndim ? ((double *)(pop->array))[0] : 0.0f;
658 v2 = pop->size>1 ? ((double *)(pop->array))[1] : v1;
659
660 /* Depending on the type of the modular warp do it. Recall that the
661 code for the warp, was stored in the 'status' element of the data
662 structure.*/
663 switch(pop->status)
664 {
665 case UI_KEY_ALIGN:
666 ui_matrix_make_align(p, module);
667 break;
668
669 case UI_KEY_ROTATE:
670 s = sin( v1 * M_PI / 180 );
671 c = cos( v1 * M_PI / 180 );
672 module[0]=c; module[1]=-1.0f*s; module[2]=0.0f;
673 module[3]=s; module[4]=c; module[5]=0.0f;
674 module[6]=0.0f; module[7]=0.0f; module[8]=1.0f;
675 break;
676
677 case UI_KEY_SCALE:
678 module[0]=v1; module[1]=0.0f; module[2]=0.0f;
679 module[3]=0.0f; module[4]=v2; module[5]=0.0f;
680 module[6]=0.0f; module[7]=0.0f; module[8]=1.0f;
681 break;
682
683 case UI_KEY_FLIP:
684 if ( v1==1.0f && v2==0.0f )
685 {
686 module[0]=1.0f; module[1]=0.0f;
687 module[3]=0.0f; module[4]=-1.0f;
688 }
689 else if ( v1==0.0f && v2==1.0f )
690 {
691 module[0]=-1.0f; module[1]=0.0f;
692 module[3]=0.0f; module[4]=1.0f;
693 }
694 else if ( v1==1.0f && v2==1.0f )
695 {
696 module[0]=-1.0f; module[1]=0.0f;
697 module[3]=0.0f; module[4]=-1.0f;
698 }
699 else /* When both are zero, just in case! */
700 {
701 module[0]=1.0f; module[1]=0.0f;
702 module[3]=0.0f; module[4]=1.0f;
703 }
704 module[2]=0.0f;
705 module[5]=0.0f;
706 module[6]=0.0f; module[7]=0.0f; module[8]=1.0f;
707 break;
708
709 case UI_KEY_SHEAR:
710 module[0]=1.0f; module[1]=v1; module[2]=0.0f;
711 module[3]=v2; module[4]=1.0f; module[5]=0.0f;
712 module[6]=0.0f; module[7]=0.0f; module[8]=1.0f;
713 break;
714
715 case UI_KEY_TRANSLATE:
716 /* Translate cannot be called with '--centeroncorner'. */
717 if(p->centeroncorner)
718 error(EXIT_FAILURE, 0, "'--translate' and '--centeroncorner' "
719 "(which is a type of translation) cannot be called "
720 "together. To achieve the effect of --centeroncorner, "
721 "start the warp steps with a translation of 0.5 to move "
722 "the coordinate center to the corner of a pixel in each "
723 "dimension");
724
725 /* Build the translation matrix. */
726 module[0]=1.0f; module[1]=0.0f; module[2]=v1;
727 module[3]=0.0f; module[4]=1.0f; module[5]=v2;
728 module[6]=0.0f; module[7]=0.0f; module[8]=1.0f;
729 break;
730
731 case UI_KEY_PROJECT:
732 module[0]=1.0f; module[1]=0.0f; module[2]=0.0f;
733 module[3]=0.0f; module[4]=1.0f; module[5]=0.0f;
734 module[6]=v1; module[7]=v2; module[8]=1.0f;
735 break;
736
737 default:
738 error(EXIT_FAILURE, 0, "a bug! the code %d is not recognized as "
739 "a valid modular warp in 'ui_matrix_from_modular', this is "
740 "not your fault, something in the programming has gone "
741 "wrong. Please contact us at %s so we can correct it",
742 pop->status, PACKAGE_BUGREPORT);
743 }
744
745 /* Multiply the main matrix with this modular matrix. */
746 ui_matrix_inplacw_multiply(p->matrix->array, module);
747
748 /* Clean up. */
749 gal_data_free(pop);
750 }
751 }
752
753
754
755
756
757 static void
ui_matrix_center_on_corner(struct warpparams * p)758 ui_matrix_center_on_corner(struct warpparams *p)
759 {
760 double *b, *d, *df;
761 double before[9]={1,0,0.5,0,1,0.5,0,0,1};
762 double after[9]={1,0,-0.5,0,1,-0.5,0,0,1};
763
764 /* Shift the matrix by +0.5 so the coordinate center lies at the bottom
765 left corner of the first pixel. Note that the updated values are
766 written into the first argument of the function.*/
767 ui_matrix_inplacw_multiply(before, p->matrix->array);
768
769 /* Translate them back into the proper FITS center. */
770 ui_matrix_inplacw_multiply(before, after);
771
772 /* The final matrix is in 'before', so put its values into the output
773 matrix. */
774 b = before;
775 df = (d=p->matrix->array) + p->matrix->size;
776 do *d=*b++; while(++d<df);
777 }
778
779
780
781
782
783 static void
ui_matrix_finalize(struct warpparams * p)784 ui_matrix_finalize(struct warpparams *p)
785 {
786 double *d, *df, *inv;
787
788 /* If a matrix string is not given, the use the modular warpings. */
789 if(p->matrix)
790 ui_matrix_prepare_raw(p);
791 else if (p->modularll)
792 ui_matrix_from_modular(p);
793 else
794 ui_error_no_warps();
795
796 /* If the user has asked for it, set the coordinate center on the corner
797 of the first pixel. */
798 if(p->centeroncorner) ui_matrix_center_on_corner(p);
799
800 /* Check if there are any non-normal numbers in the matrix: */
801 df=(d=p->matrix->array)+p->matrix->size;
802 do
803 if(!isfinite(*d++))
804 {
805 ui_matrix_print(p->matrix->array);
806 error(EXIT_FAILURE, 0, "%f is not a 'normal' number in the "
807 "input matrix shown above", *(d-1));
808 }
809 while(d<df);
810
811 /* Check if the determinant is not zero: */
812 d=p->matrix->array;
813 if( d[0]*d[4]*d[8] + d[1]*d[5]*d[6] + d[2]*d[3]*d[7]
814 - d[2]*d[4]*d[6] - d[1]*d[3]*d[8] - d[0]*d[5]*d[7] == 0 )
815 error(EXIT_FAILURE, 0, "the determinant of the given matrix "
816 "is zero");
817
818 /* Not yet implemented: Check if the transformation is spatially
819 invariant, in other words, if it differs between differet regions of
820 the output. If it doesn't we can use this information for a more
821 efficient processing. */
822
823 /* Make the inverse matrix: */
824 inv=p->inverse=gal_pointer_allocate(GAL_TYPE_FLOAT64, 9, 0, __func__,
825 "p->inverse");
826 inv[0] = d[4]*d[8] - d[5]*d[7];
827 inv[1] = d[2]*d[7] - d[1]*d[8];
828 inv[2] = d[1]*d[5] - d[2]*d[4];
829 inv[3] = d[5]*d[6] - d[3]*d[8];
830 inv[4] = d[0]*d[8] - d[2]*d[6];
831 inv[5] = d[2]*d[3] - d[0]*d[5];
832 inv[6] = d[3]*d[7] - d[4]*d[6];
833 inv[7] = d[1]*d[6] - d[0]*d[7];
834 inv[8] = d[0]*d[4] - d[1]*d[3];
835 /* Just for a test:
836 {
837 size_t i;
838 printf("\nInput matrix:");
839 for(i=0;i<9;++i) { if(i%3==0) printf("\n"); printf("%-10.5f", d[i]); }
840 printf("\n-----------\n");
841 printf("Inverse matrix:");
842 for(i=0;i<9;++i) { if(i%3==0) printf("\n"); printf("%-10.5f", inv[i]); }
843 printf("\n\n");
844 }
845 */
846 }
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867 /**************************************************************/
868 /************ General preparations ****************/
869 /**************************************************************/
870 /* When only one transformation is required, set the suffix for automatic
871 output to more meaningful string. */
872 char *
ui_set_suffix(struct warpparams * p)873 ui_set_suffix(struct warpparams *p)
874 {
875 /* A small independent sanity check: we either need a matrix or at least
876 one modular warping. */
877 if(p->matrix==NULL && p->modularll==NULL) ui_error_no_warps();
878
879 /* We only want the more meaningful suffix when the list is defined AND
880 when its only has one node (the 'next' element is NULL). */
881 if(p->matrix==NULL && p->modularll->next==NULL)
882 switch(p->modularll->status)
883 {
884 case UI_KEY_ALIGN:
885 return "_aligned.fits";
886
887 case UI_KEY_ROTATE:
888 return "_rotated.fits";
889
890 case UI_KEY_SCALE:
891 return "_scaled.fits";
892
893 case UI_KEY_FLIP:
894 return "_flipped.fits";
895
896 case UI_KEY_SHEAR:
897 return "_sheared.fits";
898
899 case UI_KEY_TRANSLATE:
900 return "_translated.fits";
901
902 case UI_KEY_PROJECT:
903 return "_projected.fits";
904
905 default:
906 error(EXIT_FAILURE, 0, "a bug! please contact us at %s so we can "
907 "fix the problem. The modular warp code %d is not recognized "
908 "in 'ui_set_suffix'", PACKAGE_BUGREPORT, p->modularll->status);
909 return NULL;
910 }
911 else
912 return "_warped.fits";
913 }
914
915
916
917
918
919 static void
ui_preparations(struct warpparams * p)920 ui_preparations(struct warpparams *p)
921 {
922 /* Set the output name. This needs to be done before 'ui_finalize_matrix'
923 because that function will free the linked list of modular warpings
924 which we will need to determine the suffix if no output name is
925 specified. */
926 if(p->cp.output)
927 gal_checkset_writable_remove(p->cp.output, 0, p->cp.dontdelete);
928 else
929 p->cp.output=gal_checkset_automatic_output(&p->cp, p->inputname,
930 ui_set_suffix(p));
931
932 /* Prepare the final warping matrix. */
933 ui_matrix_finalize(p);
934 }
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954 /**************************************************************/
955 /************ Set the parameters *************/
956 /**************************************************************/
957
958 void
ui_read_check_inputs_setup(int argc,char * argv[],struct warpparams * p)959 ui_read_check_inputs_setup(int argc, char *argv[], struct warpparams *p)
960 {
961 double *matrix;
962 struct gal_options_common_params *cp=&p->cp;
963
964
965 /* Include the parameters necessary for argp from this program ('args.h')
966 and for the common options to all Gnuastro ('commonopts.h'). We want
967 to directly put the pointers to the fields in 'p' and 'cp', so we are
968 simply including the header here to not have to use long macros in
969 those headers which make them hard to read and modify. This also helps
970 in having a clean environment: everything in those headers is only
971 available within the scope of this function. */
972 #include <gnuastro-internal/commonopts.h>
973 #include "args.h"
974
975
976 /* Initialize the options and necessary information. */
977 ui_initialize_options(p, program_options, gal_commonopts_options);
978
979
980 /* Read the command-line options and arguments. */
981 errno=0;
982 if(argp_parse(&thisargp, argc, argv, 0, 0, p))
983 error(EXIT_FAILURE, errno, "parsing arguments");
984
985
986 /* Read the configuration files and set the common values. */
987 gal_options_read_config_set(&p->cp);
988
989
990 /* Print the option values if asked. Note that this needs to be done
991 after the option checks so un-sane values are not printed in the
992 output state. */
993 gal_options_print_state(&p->cp);
994
995
996 /* Prepare all the options as FITS keywords to write in output later. */
997 gal_options_as_fits_keywords(&p->cp);
998
999
1000 /* Check that the options and arguments fit well with each other. Note
1001 that arguments don't go in a configuration file. So this test should
1002 be done after (possibly) printing the option values. */
1003 ui_check_options_and_arguments(p);
1004
1005
1006 /* Read/allocate all the necessary starting arrays. */
1007 ui_preparations(p);
1008
1009
1010 /* Everything is ready, notify the user of the program starting. */
1011 if(!p->cp.quiet)
1012 {
1013 matrix=p->matrix->array;
1014 printf(PROGRAM_NAME" "PACKAGE_VERSION" started on %s",
1015 ctime(&p->rawtime));
1016 printf(" Using %zu CPU thread%s\n", p->cp.numthreads,
1017 p->cp.numthreads==1 ? "." : "s.");
1018 printf(" Input: %s (hdu: %s)\n", p->inputname, p->cp.hdu);
1019 printf(" matrix:"
1020 "\n\t%.4f %.4f %.4f"
1021 "\n\t%.4f %.4f %.4f"
1022 "\n\t%.4f %.4f %.4f\n",
1023 matrix[0], matrix[1], matrix[2],
1024 matrix[3], matrix[4], matrix[5],
1025 matrix[6], matrix[7], matrix[8]);
1026 }
1027 }
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048 /**************************************************************/
1049 /************ Free allocated, report *************/
1050 /**************************************************************/
1051 void
ui_free_report(struct warpparams * p,struct timeval * t1)1052 ui_free_report(struct warpparams *p, struct timeval *t1)
1053 {
1054 /* Free the allocated arrays: */
1055 free(p->cp.hdu);
1056 free(p->cp.output);
1057 gal_data_free(p->input);
1058 gal_data_free(p->matrix);
1059 if(p->pixelscale) free(p->pixelscale);
1060 if(p->inwcsmatrix) free(p->inwcsmatrix);
1061
1062 /* Report how long the operation took. */
1063 if(!p->cp.quiet)
1064 gal_timing_report(t1, PROGRAM_NAME" finished in: ", 0);
1065 }
1066