1 /*********************************************************************
2 Segment - Segment initial labels based on signal structure.
3 Segment 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) 2018-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/binary.h>
35 #include <gnuastro/threads.h>
36 #include <gnuastro/dimension.h>
37 #include <gnuastro/statistics.h>
38
39 #include <gnuastro-internal/timing.h>
40 #include <gnuastro-internal/options.h>
41 #include <gnuastro-internal/checkset.h>
42 #include <gnuastro-internal/fixedstringmacros.h>
43
44 #include "main.h"
45
46 #include "ui.h"
47 #include "authors-cite.h"
48
49
50
51
52
53 /**************************************************************/
54 /********* Argp necessary global entities ************/
55 /**************************************************************/
56 /* Definition parameters for the Argp: */
57 const char *
58 argp_program_version = PROGRAM_STRING "\n"
59 GAL_STRINGS_COPYRIGHT
60 "\n\nWritten/developed by "PROGRAM_AUTHORS;
61
62 const char *
63 argp_program_bug_address = PACKAGE_BUGREPORT;
64
65 static char
66 args_doc[] = "ASTRdata";
67
68 const char
69 doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will segment an initially "
70 "labeled region based on structure with the signal. It will first find "
71 "true clumps (local maxima), estimate which ones have strong connections, "
72 "and then grow them to cover the full area of each detection.\n"
73 GAL_STRINGS_MORE_HELP_INFO
74 /* After the list of options: */
75 "\v"
76 PACKAGE_NAME" home page: "PACKAGE_URL;
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97 /**************************************************************/
98 /********* Initialize & Parse command-line **************/
99 /**************************************************************/
100 static void
ui_initialize_options(struct segmentparams * p,struct argp_option * program_options,struct argp_option * gal_commonopts_options)101 ui_initialize_options(struct segmentparams *p,
102 struct argp_option *program_options,
103 struct argp_option *gal_commonopts_options)
104 {
105 size_t i;
106 struct gal_options_common_params *cp=&p->cp;
107
108
109 /* Set the necessary common parameters structure. */
110 cp->program_struct = p;
111 cp->poptions = program_options;
112 cp->program_name = PROGRAM_NAME;
113 cp->program_exec = PROGRAM_EXEC;
114 cp->program_bibtex = PROGRAM_BIBTEX;
115 cp->program_authors = PROGRAM_AUTHORS;
116 cp->numthreads = gal_threads_number();
117 cp->coptions = gal_commonopts_options;
118
119 p->medstd = NAN;
120 p->minstd = NAN;
121 p->maxstd = NAN;
122 p->snquant = NAN;
123 p->clumpsnthresh = NAN;
124
125 /* Modify common options. */
126 for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
127 {
128 /* Select individually. */
129 switch(cp->coptions[i].key)
130 {
131 case GAL_OPTIONS_KEY_HDU:
132 cp->coptions[i].doc="HDU containing values (science image).";
133 break;
134
135 case GAL_OPTIONS_KEY_LOG:
136 case GAL_OPTIONS_KEY_TYPE:
137 case GAL_OPTIONS_KEY_SEARCHIN:
138 case GAL_OPTIONS_KEY_IGNORECASE:
139 case GAL_OPTIONS_KEY_STDINTIMEOUT:
140 cp->coptions[i].flags=OPTION_HIDDEN;
141 break;
142
143 case GAL_OPTIONS_KEY_TILESIZE:
144 case GAL_OPTIONS_KEY_MINMAPSIZE:
145 case GAL_OPTIONS_KEY_NUMCHANNELS:
146 case GAL_OPTIONS_KEY_INTERPNUMNGB:
147 case GAL_OPTIONS_KEY_REMAINDERFRAC:
148 cp->coptions[i].mandatory=GAL_OPTIONS_MANDATORY;
149 break;
150
151 case GAL_OPTIONS_KEY_TABLEFORMAT:
152 cp->coptions[i].mandatory=GAL_OPTIONS_MANDATORY;
153 cp->coptions[i].doc="'txt', 'fits-ascii', 'fits-binary'.";
154 break;
155 }
156 }
157 }
158
159
160
161
162
163 /* Parse a single option: */
164 error_t
parse_opt(int key,char * arg,struct argp_state * state)165 parse_opt(int key, char *arg, struct argp_state *state)
166 {
167 struct segmentparams *p = state->input;
168
169 /* Pass 'ygal_options_common_params' into the child parser. */
170 state->child_inputs[0] = &p->cp;
171
172 /* In case the user incorrectly uses the equal sign (for example
173 with a short format or with space in the long format, then 'arg'
174 start with (if the short version was called) or be (if the long
175 version was called with a space) the equal sign. So, here we
176 check if the first character of arg is the equal sign, then the
177 user is warned and the program is stopped: */
178 if(arg && arg[0]=='=')
179 argp_error(state, "incorrect use of the equal sign ('='). For short "
180 "options, '=' should not be used and for long options, "
181 "there should be no space between the option, equal sign "
182 "and value");
183
184 /* Set the key to this option. */
185 switch(key)
186 {
187 /* Read the non-option tokens (arguments): */
188 case ARGP_KEY_ARG:
189 /* The user may give a shell variable that is empty! In that case
190 'arg' will be an empty string! We don't want to account for such
191 cases (and give a clear error that no input has been given). */
192 if(p->inputname)
193 argp_error(state, "only one argument (input file) should be given");
194 else
195 if(arg[0]!='\0') p->inputname=arg;
196 break;
197
198 /* This is an option, set its value. */
199 default:
200 return gal_options_set_from_key(key, arg, p->cp.poptions, &p->cp);
201 }
202
203 return 0;
204 }
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225 /**************************************************************/
226 /*************** Sanity Check *******************/
227 /**************************************************************/
228 /* Read and check ONLY the options. When arguments are involved, do the
229 check in 'ui_check_options_and_arguments'. */
230 static void
ui_read_check_only_options(struct segmentparams * p)231 ui_read_check_only_options(struct segmentparams *p)
232 {
233 /* If the full area is to be used as a single detection, we can't find
234 the S/N value from the un-detected regions, so the user must have
235 given the 'clumpsnthresh' option. */
236 if( p->detectionname
237 && !strcmp(p->detectionname, DETECTION_ALL)
238 && isnan(p->clumpsnthresh) )
239 error(EXIT_FAILURE, 0, "'--clumpsnthresh' ('-%c') not given.\n\n"
240 "When '--detection=all' (the whole input dataset is assumed to "
241 "be a detection), Segment can't use the undetected pixels to find "
242 "the signal-to-noise ratio of true clumps. Therefore it is "
243 "mandatory to provide a signal-to-noise ratio manually",
244 UI_KEY_CLUMPSNTHRESH);
245
246 /* If the convolved HDU is given. */
247 if(p->convolvedname && p->chdu==NULL)
248 error(EXIT_FAILURE, 0, "no value given to '--convolvedhdu'. When the "
249 "'--convolved' option is called (to specify a convolved dataset "
250 "and avoid convolution) it is mandatory to also specify a HDU "
251 "for it");
252
253 /* For the options that make tables, the table format option is
254 mandatory. */
255 if( p->checksn && p->cp.tableformat==0 )
256 error(EXIT_FAILURE, 0, "'--tableformat' is necessary with the "
257 "'--checksn' option.\n"
258 "Please see description for '--tableformat' after running the "
259 "following command for more information (use 'SPACE' to go down "
260 "the page and 'q' to return to the command-line):\n\n"
261 " $ info gnuastro \"Input Output options\"");
262
263 /* Kernel checks. */
264 if(p->kernelname && strcmp(p->kernelname, UI_NO_CONV_KERNEL_NAME))
265 {
266 /* Check if it exists. */
267 gal_checkset_check_file(p->kernelname);
268
269 /* If its FITS, see if a HDU has been provided. */
270 if( gal_fits_file_recognized(p->kernelname) && p->khdu==NULL )
271 error(EXIT_FAILURE, 0, "no HDU specified for kernel. When the "
272 "kernel is a FITS file, a HDU must also be specified. You "
273 "can use the '--khdu' option and give it the HDU number "
274 "(starting from zero), extension name, or anything "
275 "acceptable by CFITSIO");
276 }
277
278 /* If the S/N quantile is less than 0.1 (an arbitrary small value), this
279 is probably due to forgetting that this is the purity level
280 (higher-is-better), not the contamination level
281 (lower-is-better). This actually happened in a few cases: where we
282 wanted a false detection rate of 0.0001 (a super-high value!), and
283 instead of inputing 0.9999, we mistakenly gave '--snquant' a value of
284 '0.0001'. We were thus fully confused with the output (an extremely
285 low value) and thought its a bug, while it wasn't! */
286 if(p->snquant<0.1)
287 fprintf(stderr, "\nWARNING: Value of '--snquant' ('-c') is %g. Note "
288 "that this is not a contamination rate (where lower is "
289 "better), it is a purity rate (where higher is better). If you "
290 "intentionally asked for such a low purity level, please "
291 "ignore this warning\n\n", p->snquant);
292 }
293
294
295
296
297
298 static void
ui_check_options_and_arguments(struct segmentparams * p)299 ui_check_options_and_arguments(struct segmentparams *p)
300 {
301 /* Make sure an input file name was given and if it was a FITS file, that
302 a HDU is also given. */
303 if(p->inputname)
304 {
305 /* Check if it exists. */
306 gal_checkset_check_file(p->inputname);
307
308 /* If it is FITS, a HDU is also mandatory. */
309 if( gal_fits_file_recognized(p->inputname) && p->cp.hdu==NULL )
310 error(EXIT_FAILURE, 0, "no HDU specified. When the input is a FITS "
311 "file, a HDU must also be specified, you can use the '--hdu' "
312 "('-h') option and give it the HDU number (starting from "
313 "zero), extension name, or anything acceptable by CFITSIO");
314
315 }
316 else
317 error(EXIT_FAILURE, 0, "no input file is specified");
318 }
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339 /**************************************************************/
340 /*************** Preparations *******************/
341 /**************************************************************/
342 static void
ui_set_used_names(struct segmentparams * p)343 ui_set_used_names(struct segmentparams *p)
344 {
345 p->useddetectionname = p->detectionname ? p->detectionname : p->inputname;
346
347 p->usedstdname = ( p->stdname
348 ? p->stdname
349 : ( ( p->detectionname
350 && strcmp(p->detectionname, DETECTION_ALL) )
351 ? p->detectionname
352 : p->inputname ) );
353 }
354
355
356
357
358
359 static void
ui_set_output_names(struct segmentparams * p)360 ui_set_output_names(struct segmentparams *p)
361 {
362 char *output=p->cp.output;
363 char *basename = output ? output : p->inputname;
364
365 /* Main program output. */
366 if(output)
367 {
368 /* Delete the file if it already exists. */
369 gal_checkset_writable_remove(p->cp.output, 0, p->cp.dontdelete);
370
371 /* When the output name is given (possibly with directory
372 information), the check images will also be put in that same
373 directory. */
374 p->cp.keepinputdir=1;
375 }
376 else
377 p->cp.output=gal_checkset_automatic_output(&p->cp, p->inputname,
378 "_segmented.fits");
379
380 /* Tile check. */
381 if(p->cp.tl.checktiles)
382 p->cp.tl.tilecheckname=gal_checkset_automatic_output(&p->cp, basename,
383 "_tiles.fits");
384
385 /* Clump S/N values. */
386 if(p->checksn)
387 {
388 p->clumpsn_s_name=gal_checkset_automatic_output(&p->cp, basename,
389 ( p->cp.tableformat==GAL_TABLE_FORMAT_TXT
390 ? "_clumpsn_sky.txt" : "_clumpsn.fits") );
391 p->clumpsn_d_name=gal_checkset_automatic_output(&p->cp, basename,
392 ( p->cp.tableformat==GAL_TABLE_FORMAT_TXT
393 ? "_clumpsn_det.txt" : "_clumpsn.fits") );
394 }
395
396 /* Segmentation steps. */
397 if(p->checksegmentation)
398 p->segmentationname=gal_checkset_automatic_output(&p->cp, basename,
399 "_segcheck.fits");
400 }
401
402
403
404
405
406 static void
ui_prepare_inputs(struct segmentparams * p)407 ui_prepare_inputs(struct segmentparams *p)
408 {
409 int32_t *i, *ii;
410 gal_data_t *maxd, *ccin, *blankflag, *ccout=NULL;
411
412 /* Read the input as a single precision floating point dataset. */
413 p->input = gal_array_read_one_ch_to_type(p->inputname, p->cp.hdu,
414 NULL, GAL_TYPE_FLOAT32,
415 p->cp.minmapsize,
416 p->cp.quietmmap);
417 p->input->wcs = gal_wcs_read(p->inputname, p->cp.hdu,
418 p->cp.wcslinearmatrix, 0, 0,
419 &p->input->nwcs);
420 p->input->ndim=gal_dimension_remove_extra(p->input->ndim,
421 p->input->dsize,
422 p->input->wcs);
423
424 /* Set the name. */
425 if(p->input->name) free(p->input->name);
426 gal_checkset_allocate_copy("INPUT", &p->input->name);
427
428
429 /* Check for blank values to help later processing. */
430 gal_blank_present(p->input, 1);
431
432
433 /* Segment currently only works on 2D datasets (images). */
434 if(p->input->ndim!=2 && p->input->ndim!=3)
435 error(EXIT_FAILURE, 0, "%s (hdu: %s) has %zu dimensions but Segment "
436 "can only operate on 2D (images) or 3D (cube) datasets",
437 p->inputname, p->cp.hdu, p->input->ndim);
438
439
440 /* If a convolved image is given, read it. */
441 if(p->convolvedname)
442 {
443 /* Read the input convolved image. */
444 p->conv = gal_array_read_one_ch_to_type(p->convolvedname, p->chdu,
445 NULL, GAL_TYPE_FLOAT32,
446 p->cp.minmapsize,
447 p->cp.quietmmap);
448 p->conv->ndim=gal_dimension_remove_extra(p->conv->ndim,
449 p->conv->dsize,
450 p->conv->wcs);
451 p->conv->wcs=gal_wcs_copy(p->input->wcs);
452
453 /* Make sure it is the same size as the input. */
454 if( gal_dimension_is_different(p->input, p->conv) )
455 error(EXIT_FAILURE, 0, "%s (hdu %s), given to '--convolved' and "
456 "'--chdu', is not the same size as the input (%s, hdu: %s)",
457 p->convolvedname, p->chdu, p->inputname, p->cp.hdu);
458 }
459
460
461 /* Read the detected label image and check its size. When the user gives
462 '--detection=all', then the whole input is assumed to be a single
463 detection. */
464 if( strcmp(p->useddetectionname, DETECTION_ALL) )
465 {
466 /* Read the dataset into memory. */
467 p->olabel = gal_array_read_one_ch(p->useddetectionname, p->dhdu,
468 NULL, p->cp.minmapsize,
469 p->cp.quietmmap);
470 p->olabel->ndim=gal_dimension_remove_extra(p->olabel->ndim,
471 p->olabel->dsize, NULL);
472 if( gal_dimension_is_different(p->input, p->olabel) )
473 error(EXIT_FAILURE, 0, "'%s' (hdu: %s) and '%s' (hdu: %s) have a"
474 "different dimension/size", p->useddetectionname, p->dhdu,
475 p->inputname, p->cp.hdu);
476
477 /* Make sure the detected labels are not floating point. */
478 if(p->olabel->type==GAL_TYPE_FLOAT32
479 || p->olabel->type==GAL_TYPE_FLOAT64)
480 error(EXIT_FAILURE, 0, "%s (hdu: %s) has a '%s' type. The detection "
481 "(labeled) map must have an integer type (labels/classes can "
482 "only be integers). If the pixel values are integers, but only "
483 "the numerical type of the image is floating-point, you can "
484 "use the command below to convert it to a 32-bit (signed) "
485 "integer type:\n\n"
486 " $ astarithmetic %s int32 -h%s\n\n", p->useddetectionname,
487 p->dhdu, gal_type_name(p->olabel->type, 1),
488 p->useddetectionname, p->dhdu);
489
490 /* If the input has blank values, set them to blank values in the
491 labeled image too. It doesn't matter if the labeled image has
492 blank pixels that aren't blank on the input image. */
493 if(gal_blank_present(p->input, 1))
494 {
495 blankflag=gal_blank_flag(p->input);
496 gal_blank_flag_apply(p->olabel, blankflag);
497 gal_data_free(blankflag);
498 }
499
500 /* Get the maximum value of the input (total number of labels if they
501 are separate). If the maximum is 1 (the image is a binary image),
502 then apply the connected components algorithm to separate the
503 connected regions. The user is allowed to supply a simple binary
504 image.*/
505 maxd=gal_statistics_maximum(p->olabel);
506 maxd=gal_data_copy_to_new_type_free(maxd, GAL_TYPE_INT64);
507 p->numdetections = *((uint64_t *)(maxd->array));
508 if( p->numdetections == 1 )
509 {
510 ccin=gal_data_copy_to_new_type_free(p->olabel, GAL_TYPE_UINT8);
511 p->numdetections=gal_binary_connected_components(ccin, &ccout,
512 ccin->ndim);
513 gal_data_free(ccin);
514 p->olabel=ccout;
515 }
516 else
517 p->olabel = gal_data_copy_to_new_type_free(p->olabel, GAL_TYPE_INT32);
518
519
520 /* Write the WCS into the objects dataset too. */
521 p->olabel->wcs=gal_wcs_copy(p->input->wcs);
522 }
523 else
524 {
525 /* Set the total number of detections to 1. */
526 p->numdetections=1;
527
528 /* Allocate the array. */
529 p->olabel=gal_data_alloc(NULL, GAL_TYPE_INT32, p->input->ndim,
530 p->input->dsize, p->input->wcs, 0,
531 p->cp.minmapsize, p->cp.quietmmap,
532 NULL, NULL, NULL);
533
534 /* Initialize it to 1. */
535 ii=(i=p->olabel->array)+p->olabel->size; do *i++=1; while(i<ii);
536 }
537 }
538
539
540
541
542
543 /* Prepare the kernel, either from a file, or from the default arrays
544 available in the headers. The default kernels were created as
545 follows. */
546 static void
ui_prepare_kernel(struct segmentparams * p)547 ui_prepare_kernel(struct segmentparams *p)
548 {
549 float *f, *ff, *k;
550 size_t ndim=p->input->ndim;
551
552 /* Import the default kernel. */
553 #include "kernel-2d.h"
554 #include "kernel-3d.h"
555
556 /* If a kernel file is given, then use it. Otherwise, use the default
557 kernel. */
558 if(p->kernelname)
559 {
560 /* Read the kernel. */
561 if( strcmp(p->kernelname, UI_NO_CONV_KERNEL_NAME) )
562 {
563 /* Read the kernel into memory. */
564 p->kernel=gal_fits_img_read_kernel(p->kernelname, p->khdu,
565 p->cp.minmapsize,
566 p->cp.quietmmap);
567 p->kernel->ndim=gal_dimension_remove_extra(p->kernel->ndim,
568 p->kernel->dsize,
569 NULL);
570
571 /* Make sure it has the same dimensions as the input. */
572 if( p->kernel->ndim != p->input->ndim )
573 error(EXIT_FAILURE, 0, "%s (hdu %s): is %zuD, however, %s (%s) "
574 "is a %zuD dataset", p->kernelname, p->khdu,
575 p->kernel->ndim, p->inputname, p->cp.hdu, p->input->ndim);
576 }
577 else
578 p->kernel=NULL;
579 }
580 else
581 {
582 /* Allocate space for the kernel (we don't want to use the statically
583 allocated array. */
584 p->kernel=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, p->input->ndim,
585 ndim==2 ? kernel_2d_dsize : kernel_3d_dsize,
586 NULL, 0, p->cp.minmapsize, p->cp.quietmmap,
587 NULL, NULL, NULL);
588
589 /* Copy the staticly allocated default array into `p->kernel'. */
590 k = ndim==2 ? kernel_2d : kernel_3d;
591 ff = (f=p->kernel->array) + p->kernel->size;
592 do *f=*k++; while(++f<ff);
593 }
594 }
595
596
597
598
599
600 /* Set up the tessellation. */
601 static void
ui_prepare_tiles(struct segmentparams * p)602 ui_prepare_tiles(struct segmentparams *p)
603 {
604 gal_data_t *check;
605 struct gal_tile_two_layer_params *tl=&p->cp.tl, *ltl=&p->ltl;
606
607
608 /* Check the tile parameters for the small tile sizes and make the tile
609 structure. */
610 gal_tile_full_sanity_check(p->inputname, p->cp.hdu, p->input, tl);
611 gal_tile_full_two_layers(p->input, tl);
612 gal_tile_full_permutation(tl);
613
614
615 /* Make the large tessellation, except for the size, the rest of the
616 parameters are the same as the small tile sizes. */
617 ltl->numchannels = tl->numchannels;
618 ltl->remainderfrac = tl->remainderfrac;
619 ltl->workoverch = tl->workoverch;
620 ltl->checktiles = tl->checktiles;
621 ltl->oneelempertile = tl->oneelempertile;
622 gal_tile_full_sanity_check(p->inputname, p->cp.hdu, p->input, ltl);
623 gal_tile_full_two_layers(p->input, ltl);
624 gal_tile_full_permutation(ltl);
625
626
627 /* If the input has blank elements, then set the appropriate flag for
628 each tile.*/
629 if( p->input->flag & GAL_DATA_FLAG_HASBLANK )
630 {
631 gal_tile_block_blank_flag(tl->tiles, p->cp.numthreads);
632 gal_tile_block_blank_flag(ltl->tiles, p->cp.numthreads);
633 }
634
635
636 /* Make the tile check image if requested. */
637 if(tl->checktiles)
638 {
639 /* Large tiles. */
640 check=gal_tile_block_check_tiles(ltl->tiles);
641 gal_fits_img_write(check, tl->tilecheckname, NULL, PROGRAM_NAME);
642 gal_data_free(check);
643
644 /* Small tiles. */
645 check=gal_tile_block_check_tiles(tl->tiles);
646 gal_fits_img_write(check, tl->tilecheckname, NULL, PROGRAM_NAME);
647 gal_data_free(check);
648
649 /* If 'continueaftercheck' hasn't been called, abort NoiseChisel. */
650 if(!p->continueaftercheck)
651 ui_abort_after_check(p, tl->tilecheckname, NULL,
652 "showing all tiles over the image");
653
654 /* Free the name. */
655 free(tl->tilecheckname);
656 tl->tilecheckname=NULL;
657 }
658 }
659
660
661
662
663
664 static void
ui_check_size(gal_data_t * base,gal_data_t * comp,size_t numtiles,char * bname,char * bhdu,char * cname,char * chdu)665 ui_check_size(gal_data_t *base, gal_data_t *comp, size_t numtiles,
666 char *bname, char *bhdu, char *cname, char *chdu)
667 {
668 if( gal_dimension_is_different(base, comp) && numtiles!=comp->size )
669 error(EXIT_FAILURE, 0, "%s (hdu: %s): doesn't have the right size "
670 "(%zu elements or pixels).\n\n"
671 "It must either be the same size as '%s' (hdu: '%s'), or "
672 "it must have the same number of elements as the total "
673 "number of tiles in the tessellation (%zu). In the latter "
674 "case, each pixel is assumed to be a fixed value for a "
675 "complete tile.\n\n"
676 "Run with '-P' to see the (tessellation) options/settings "
677 "and their values). For more information on tessellation in "
678 "Gnuastro, please run the following command (use the arrow "
679 "keys for up and down and press 'q' to return to the "
680 "command-line):\n\n"
681 " $ info gnuastro tessellation",
682 cname, chdu, comp->size, bname, bhdu, numtiles);
683 }
684
685
686
687
688
689 /* Subtract 'sky' from the input dataset depending on its size (it may be
690 the whole array or a tile-values array).. */
691 static void
ui_subtract_sky(gal_data_t * in,gal_data_t * sky,struct gal_tile_two_layer_params * tl)692 ui_subtract_sky(gal_data_t *in, gal_data_t *sky,
693 struct gal_tile_two_layer_params *tl)
694 {
695 size_t tid;
696 gal_data_t *tile;
697 float *s, *f, *ff, *skyarr=sky->array;
698
699 /* It is the same size as the input or a single value. */
700 if( gal_dimension_is_different(in, sky)==0 || sky->size==1)
701 {
702 s=sky->array;
703 ff=(f=in->array)+in->size;
704 if(sky->size==1) { if(*s!=0.0) do *f-=*s; while(++f<ff); }
705 else do *f-=*s++; while(++f<ff);
706 }
707
708 /* It is the same size as the number of tiles. */
709 else if( tl->tottiles==sky->size )
710 {
711 /* Go over all the tiles. */
712 for(tid=0; tid<tl->tottiles; ++tid)
713 {
714 /* For easy reading. */
715 tile=&tl->tiles[tid];
716
717 /* Subtract the Sky value from the input image. */
718 GAL_TILE_PARSE_OPERATE(tile, NULL, 0, 0, {*i-=skyarr[tid];});
719 }
720 }
721
722 /* The size must have been checked before, so if control reaches here, we
723 have a bug! */
724 else
725 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
726 "the problem. For some reason, the size doesn't match", __func__,
727 PACKAGE_BUGREPORT);
728 }
729
730
731
732
733
734 /* The Sky and Sky standard deviation images can be a 'oneelempertile'
735 image (only one element/pixel for a tile). So we need to do some extra
736 checks on them (after reading the tessellation). */
737 static float
ui_read_std_and_sky(struct segmentparams * p)738 ui_read_std_and_sky(struct segmentparams *p)
739 {
740 size_t one=1;
741 char *tailptr;
742 float tmpval, skyval=NAN;
743 struct gal_tile_two_layer_params *tl=&p->cp.tl;
744 gal_data_t *sky, *keys=gal_data_array_calloc(3);
745
746 /* See if the name used for the standard deviation is a filename or a
747 value. When the string is only a number (and nothing else), 'tailptr'
748 will point to the end of the string ('\0'). When the string doesn't
749 start with a number, it will point to the start of the
750 string. However, file names might also be things like '1_std.fits'. In
751 such cases, 'strtod' will return '1.0' and 'tailptr' will be
752 '_std.fits'. Thus the most robust test is to see if 'tailptr' is the
753 NULL string character. */
754 tmpval=strtod(p->usedstdname, &tailptr);
755 if(*tailptr=='\0')
756 {
757 /* Allocate the dataset to keep the value and write it in. */
758 p->std=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0,
759 -1, 1, NULL, NULL, NULL);
760 *(float *)(p->std->array) = tmpval;
761 }
762 else
763 {
764 /* Make sure a HDU is also given. */
765 if(p->stdhdu==NULL)
766 error(EXIT_FAILURE, 0, "no value given to '--stdhdu'.\n\n"
767 "When the Sky standard deviation is a dataset, it is mandatory "
768 "specify which HDU/extension it is present in. The file can "
769 "be specified explicitly with '--std'. If not, segment will "
770 "use the file given to '--detection'. If that is also not "
771 "called, it will look into the main input file (with no "
772 "option)");
773
774 /* Read the STD image. */
775 p->std=gal_array_read_one_ch_to_type(p->usedstdname, p->stdhdu,
776 NULL, GAL_TYPE_FLOAT32,
777 p->cp.minmapsize, p->cp.quietmmap);
778 p->std->ndim=gal_dimension_remove_extra(p->std->ndim,
779 p->std->dsize, NULL);
780
781 /* Make sure it has the correct size. */
782 ui_check_size(p->input, p->std, tl->tottiles, p->inputname, p->cp.hdu,
783 p->usedstdname, p->stdhdu);
784 }
785
786 /* When the Standard deviation dataset (not single value) is made by
787 NoiseChisel, it puts three basic statistics of the pre-interpolation
788 distribution of standard deviations in 'MEDSTD', 'MINSTD' and
789 'MAXSTD'. The 'MEDSTD' in particular is most important because it
790 can't be inferred after the interpolations and it can be useful in
791 MakeCatalog later to give a more accurate estimate of the noise
792 level. So if they are present, we will read them here and write them
793 to the STD output (which is created when '--rawoutput' is not
794 given). */
795 if(!p->rawoutput && p->std->size>1)
796 {
797 keys[0].next=&keys[1];
798 keys[1].next=&keys[2];
799 keys[2].next=NULL;
800 keys[0].array=&p->medstd; keys[0].name="MEDSTD";
801 keys[1].array=&p->minstd; keys[1].name="MINSTD";
802 keys[2].array=&p->maxstd; keys[2].name="MAXSTD";
803 keys[0].type=keys[1].type=keys[2].type=GAL_TYPE_FLOAT32;
804 gal_fits_key_read(p->usedstdname, p->stdhdu, keys, 0, 0);
805 if(keys[0].status) p->medstd=NAN;
806 if(keys[1].status) p->minstd=NAN;
807 if(keys[2].status) p->maxstd=NAN;
808 keys[0].name=keys[1].name=keys[2].name=NULL;
809 keys[0].array=keys[1].array=keys[2].array=NULL;
810 gal_data_array_free(keys, 3, 1);
811 }
812
813 /* Similar to '--std' above. */
814 if(p->skyname)
815 {
816 tmpval=strtod(p->skyname, &tailptr);
817 if(*tailptr=='\0')
818 {
819 sky=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0, -1,
820 1, NULL, NULL, NULL);
821 *(float *)(sky->array) = skyval = tmpval;
822 }
823 else
824 {
825 /* Make sure a HDU is also given. */
826 if(p->skyhdu==NULL)
827 error(EXIT_FAILURE, 0, "no value given to '--skyhdu'.\n\n"
828 "When the Sky is a dataset, it is mandatory specify "
829 "which HDU/extension it is present in. The file can be "
830 "specified explicitly with '--sky'. If it is a single "
831 "value, you can just pass the value to '--sky' and no "
832 "HDU will be necessary");
833
834 /* Read the Sky dataset. */
835 sky=gal_array_read_one_ch_to_type(p->skyname, p->skyhdu,
836 NULL, GAL_TYPE_FLOAT32,
837 p->cp.minmapsize, p->cp.quietmmap);
838 sky->ndim=gal_dimension_remove_extra(sky->ndim, sky->dsize,
839 NULL);
840
841 /* Check its size. */
842 ui_check_size(p->input, sky, tl->tottiles, p->inputname, p->cp.hdu,
843 p->skyname, p->skyhdu);
844 }
845
846 /* Subtract the sky from the input. */
847 ui_subtract_sky(p->input, sky, tl);
848
849 /* If a convolved image is given, subtract the Sky from that too. */
850 if(p->conv)
851 ui_subtract_sky(p->conv, sky, tl);
852
853 /* Clean up. */
854 gal_data_free(sky);
855 }
856
857 /* Return the sky value (possibly necessary in verbose mode). */
858 return skyval;
859 }
860
861
862
863
864
865 static float
ui_preparations(struct segmentparams * p)866 ui_preparations(struct segmentparams *p)
867 {
868 /* Set the input names. */
869 ui_set_used_names(p);
870
871 /* Prepare the names of the outputs. */
872 ui_set_output_names(p);
873
874 /* Read the input datasets. */
875 ui_prepare_inputs(p);
876
877 /* If a convolved image was given, read it in. Otherwise, read the given
878 kernel. */
879 if(p->conv==NULL)
880 ui_prepare_kernel(p);
881
882 /* Prepare the tessellation. */
883 ui_prepare_tiles(p);
884
885 /* Prepare the (optional Sky, and) Sky Standard deviation image. */
886 return ui_read_std_and_sky(p);
887 }
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907 /**************************************************************/
908 /************ Set the parameters *************/
909 /**************************************************************/
910 void
ui_read_check_inputs_setup(int argc,char * argv[],struct segmentparams * p)911 ui_read_check_inputs_setup(int argc, char *argv[], struct segmentparams *p)
912 {
913 float sky;
914 char *stdunit;
915 struct gal_options_common_params *cp=&p->cp;
916
917
918 /* Include the parameters necessary for argp from this program ('args.h')
919 and for the common options to all Gnuastro ('commonopts.h'). We want
920 to directly put the pointers to the fields in 'p' and 'cp', so we are
921 simply including the header here to not have to use long macros in
922 those headers which make them hard to read and modify. This also helps
923 in having a clean environment: everything in those headers is only
924 available within the scope of this function. */
925 #include <gnuastro-internal/commonopts.h>
926 #include "args.h"
927
928
929 /* Initialize the options and necessary information. */
930 ui_initialize_options(p, program_options, gal_commonopts_options);
931
932
933 /* Read the command-line options and arguments. */
934 errno=0;
935 if(argp_parse(&thisargp, argc, argv, 0, 0, p))
936 error(EXIT_FAILURE, errno, "parsing arguments");
937
938
939 /* Read the configuration files and set the common values. */
940 gal_options_read_config_set(&p->cp);
941
942
943 /* Read the options into the program's structure, and check them and
944 their relations prior to printing. */
945 ui_read_check_only_options(p);
946
947
948 /* Print the option values if asked. Note that this needs to be done
949 after the option checks so un-sane values are not printed in the
950 output state. */
951 gal_options_print_state(&p->cp);
952
953
954 /* Prepare all the options as FITS keywords to write in output later. */
955 gal_options_as_fits_keywords(&p->cp);
956
957
958 /* Check that the options and arguments fit well with each other. Note
959 that arguments don't go in a configuration file. So this test should
960 be done after (possibly) printing the option values. */
961 ui_check_options_and_arguments(p);
962
963
964 /* Read/allocate all the necessary starting arrays. */
965 sky=ui_preparations(p);
966
967
968 /* Let the user know that processing has started. */
969 if(!p->cp.quiet)
970 {
971 /* Basic inputs. */
972 printf(PROGRAM_NAME" "PACKAGE_VERSION" started on %s",
973 ctime(&p->rawtime));
974 printf(" - Using %zu CPU thread%s\n", p->cp.numthreads,
975 p->cp.numthreads==1 ? "." : "s.");
976 printf(" - Input: %s (hdu: %s)\n", p->inputname, p->cp.hdu);
977
978 /* Sky value information. */
979 if(p->skyname)
980 {
981 if( isnan(sky) )
982 printf(" - Sky: %s (hdu: %s)\n", p->skyname, p->skyhdu);
983 else
984 printf(" - Sky: %g\n", sky);
985 }
986
987 /* Sky Standard deviation information. */
988 stdunit = p->variance ? "VAR" : "STD";
989 if(p->std->size>1)
990 printf(" - Sky %s: %s (hdu: %s)\n", stdunit, p->usedstdname,
991 p->stdhdu);
992 else
993 printf(" - Sky %s: %g\n", stdunit, *(float *)(p->std->array) );
994
995 /* Convolution information. */
996 if(p->convolvedname)
997 printf(" - Convolved input: %s (hdu: %s)\n", p->convolvedname,
998 p->chdu);
999 else
1000 {
1001 if(p->kernelname)
1002 {
1003 if( strcmp(p->kernelname, UI_NO_CONV_KERNEL_NAME) )
1004 printf(" - Kernel: %s (hdu: %s)\n", p->kernelname, p->khdu);
1005 else
1006 printf(" - No convolution requested.\n");
1007 }
1008 else
1009 printf(" - Kernel: FWHM=1.5 pixel Gaussian.\n");
1010 }
1011 printf(" - Detection: %s (hdu: %s)\n", p->useddetectionname, p->dhdu);
1012 }
1013 }
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034 /**************************************************************/
1035 /************ Free allocated, report *************/
1036 /**************************************************************/
1037 void
ui_abort_after_check(struct segmentparams * p,char * filename,char * file2name,char * description)1038 ui_abort_after_check(struct segmentparams *p, char *filename,
1039 char *file2name, char *description)
1040 {
1041 char *name;
1042
1043 if(file2name)
1044 {
1045 if( asprintf(&name, "'%s' and '%s'", filename, file2name)<0 )
1046 error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
1047 }
1048 else
1049 {
1050 if( asprintf(&name, "'%s'", filename)<0 )
1051 error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
1052 }
1053
1054 /* Let the user know that NoiseChisel is aborting. */
1055 fprintf(stderr,
1056 "------------------------------------------------\n"
1057 "%s aborted for a check\n"
1058 "------------------------------------------------\n"
1059 "%s (%s) has been created.\n\n"
1060 "If you want %s to continue its processing AND save any "
1061 "requested check outputs, please run it again with "
1062 "'--continueaftercheck'.\n"
1063 "------------------------------------------------\n",
1064 PROGRAM_NAME, name, description, PROGRAM_NAME);
1065
1066 /* Clean up. */
1067 free(name);
1068 ui_free_report(p, NULL);
1069
1070 /* Abort. */
1071 exit(EXIT_SUCCESS);
1072 }
1073
1074
1075
1076
1077
1078 void
ui_free_report(struct segmentparams * p,struct timeval * t1)1079 ui_free_report(struct segmentparams *p, struct timeval *t1)
1080 {
1081 /* Free the allocated arrays: */
1082 free(p->cp.hdu);
1083 free(p->cp.output);
1084 gal_data_free(p->std);
1085 gal_data_free(p->input);
1086 gal_data_free(p->kernel);
1087 gal_data_free(p->binary);
1088 gal_data_free(p->olabel);
1089 gal_data_free(p->clabel);
1090 if(p->khdu) free(p->khdu);
1091 if(p->chdu) free(p->chdu);
1092 if(p->dhdu) free(p->dhdu);
1093 if(p->skyhdu) free(p->skyhdu);
1094 if(p->stdhdu) free(p->stdhdu);
1095 if(p->stdname) free(p->stdname);
1096 if(p->kernelname) free(p->kernelname);
1097 if(p->detectionname) free(p->detectionname);
1098 if(p->convolvedname) free(p->convolvedname);
1099 if(p->conv!=p->input) gal_data_free(p->conv);
1100 if(p->clumpsn_s_name) free(p->clumpsn_s_name);
1101 if(p->clumpsn_d_name) free(p->clumpsn_d_name);
1102 if(p->segmentationname) free(p->segmentationname);
1103
1104 /* Print the final message. */
1105 if(!p->cp.quiet && t1)
1106 gal_timing_report(t1, PROGRAM_NAME" finished in: ", 0);
1107 }
1108