[gnuastro-commits] master cc5ab9df 1/6: dither-simulate: new script to simulate a dither pattern stack

[Mohammad Akhlaghi]

branch: master
commit cc5ab9df741e1732d76eaeefce481b928fc389d7
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    dither-simulate: new script to simulate a dither pattern stack
    
    Until now, there was no easy, high-level way for Gnuastro users to simulate
    a certain dither pattern. But this is often necessary when we are planning
    an observation.
    
    With this commit, a new script has been added to do the job. This is its
    raw format and it will be modified and documented in the next series
    commits before this branch gets merged with the 'master' branch.
---
 bin/script/Makefile.am        |   6 +
 bin/script/dither-simulate.in | 339 ++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 345 insertions(+)

diff --git a/bin/script/Makefile.am b/bin/script/Makefile.am
index 6c12f854..21791a6a 100644
--- a/bin/script/Makefile.am
+++ b/bin/script/Makefile.am
@@ -45,6 +45,7 @@ bin_SCRIPTS = astscript-fits-view \
               astscript-psf-subtract \
               astscript-sort-by-night \
               astscript-radial-profile \
+              astscript-dither-simulate \
               astscript-psf-scale-factor \
               astscript-psf-select-stars
 
@@ -57,6 +58,7 @@ EXTRA_DIST = fits-view.in \
              psf-subtract.in \
              sort-by-night.in \
              radial-profile.in \
+             dither-simulate.in \
              psf-select-stars.in \
              psf-scale-factor.in \
              fits-view.desktop.in
@@ -79,6 +81,10 @@ do_subst = sed -e 's,[@]PREFIX[@],$(exec_prefix),g' \
 
 
 ## Rules to install the scripts.
+astscript-dither-simulate: dither-simulate.in Makefile
+       $(do_subst) < $(srcdir)/dither-simulate.in > $@
+       chmod +x $@
+
 astscript-ds9-region: ds9-region.in Makefile
        $(do_subst) < $(srcdir)/ds9-region.in > $@
        chmod +x $@
diff --git a/bin/script/dither-simulate.in b/bin/script/dither-simulate.in
new file mode 100644
index 00000000..a18afec0
--- /dev/null
+++ b/bin/script/dither-simulate.in
@@ -0,0 +1,339 @@
+#!/bin/sh
+# Script to simuate a dither pattern based on an existing image
+# (accounting for its distortion).
+#
+# Copyright (C) 2023-2023 Mohammad Akhlaghi <mohammad@akhlaghi.org>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+# Exit the script in the case of failure
+set -e
+
+# 'LC_NUMERIC' is responsible for formatting numbers printed by the OS.  It
+# prevents floating points like '23,45' instead of '23.45'.
+export LC_NUMERIC=C
+
+# Default values.
+img=""
+cat=""
+width=""
+imghdu=1
+center=""
+keeptmp=0
+output=""
+racol="RA"
+deccol="DEC"
+widthinpix=0
+version=@VERSION@
+scriptname=@SCRIPT_NAME@
+output=dither-check.fits
+
+
+
+
+
+# Output of '--help'
+print_help() {
+   cat <<EOF
+Usage: $scriptname positions-cat.fits --center=1.23,4.56 \
+                   --width=2,2 --img=image.fits [OPTIONS]
+
+Given a set of dithering positions ('positions-cat.fits' in the example
+above), and an image ('image.fits' which can be from any part of the sky,
+only its distortion and orientation are important), build the exposure map
+of the output stack after applying that dither (where each pixel contains
+the number of exposures that were used in it). The field of the final stack
+can be set with the '--center' and '--width' options.
+
+$scriptname options:
+ Input:
+  -h, --hdu=STR/INT       HDU name or number of the positions catalog.
+  -i, --img=STR           Name of reference image.
+  -H, --imghdu=STR/INT    HDU name or number of the input image.
+  -r, --racol=STR/INT     Name/number of column containing RA.
+  -d, --deccol=STR/INT    Name/number of column containing Declination.
+  -C, --center=FLT,FLT    Center of output stack (in RA,Dec).
+  -w, --width=FLT,FLT     Width of output stack (in WCS).
+      --widthinpix        Interpret '--width' values in pixels.
+
+ Output:
+  -o, --output=STR        Name of finally stacked image.
+  -t, --tmpdir=STR        Directory to keep temporary files.
+  -k, --keeptmp           Keep temporal/auxiliar files.
+
+ Operating mode:
+  -h, --help              Print this help.
+      --cite              BibTeX citation for this program.
+  -q, --quiet             Don't print any extra information in stdout.
+  -V, --version           Print program version.
+
+Mandatory or optional arguments to long options are also mandatory or optional
+for any corresponfing short options.
+
+GNU Astronomy Utilities home page: http://www.gnu.org/software/gnuastro/
+
+Report bugs to bug-gnuastro@gnu.org
+EOF
+}
+
+print_version() {
+     cat <<EOF
+$scriptname (GNU Astronomy Utilities) $version
+Copyright (C) 2020-2023 Free Software Foundation, Inc.
+License GPLv3+: GNU General public license version 3 or later.
+This is free software: you are free to change and redistribute it.
+There is NO WARRANTY, to the extent permitted by law.
+
+Written/developed by Mohammad Akhlaghi, Raul Infante-Sainz and Sepideh 
Eskandarlou.
+EOF
+}
+
+print_citation() {
+    empty="" # needed for the ascii art!
+    cat <<EOF
+
+Thank you for using $scriptname (GNU Astronomy Utilities) $version
+
+Citations and acknowledgement are vital for the continued work on Gnuastro.
+
+Please cite the following record(s) and add the acknowledgement statement 
below in your work to support us. Please note that different Gnuastro programs 
may have different corresponding papers. Hence, please check all the programs 
you used. Don't forget to also include the version as shown above for 
reproducibility.
+
+First paper introducing Gnuastro
+--------------------------------
+  @ARTICLE{gnuastro,
+     author = {{Akhlaghi}, M. and {Ichikawa}, T.},
+      title = "{Noise-based Detection and Segmentation of Nebulous Objects}",
+    journal = {ApJS},
+  archivePrefix = "arXiv",
+     eprint = {1505.01664},
+   primaryClass = "astro-ph.IM",
+       year = 2015,
+      month = sep,
+     volume = 220,
+        eid = {1},
+      pages = {1},
+        doi = {10.1088/0067-0049/220/1/1},
+     adsurl = {https://ui.adsabs.harvard.edu/abs/2015ApJS..220....1A},
+    adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+  }
+
+Acknowledgement
+---------------
+This work was partly done using GNU Astronomy Utilities (Gnuastro, 
ascl.net/1801.009) version $version. Work on Gnuastro has been funded by the 
Japanese Ministry of Education, Culture, Sports, Science, and Technology (MEXT) 
scholarship and its Grant-in-Aid for Scientific Research (21244012, 24253003), 
the European Research Council (ERC) advanced grant 339659-MUSICOS, the Spanish 
Ministry of Economy and Competitiveness (MINECO, grant number AYA2016-76219-P) 
and the NextGenerationEU grant  [...]
+                                               ,
+                                              {|'--.
+                                             {{\    \ $empty
+      Many thanks from all                   |/\`'--./=.
+      Gnuastro developers!                   \`\.---' \`\\
+                                                  |\  ||
+                                                  | |//
+                                                   \//_/|
+                                                   //\__/
+                                                  //
+                   (http://www.chris.com/ascii/) |/
+
+EOF
+}
+
+on_off_option_error() {
+    if [ x"$2" = x ]; then
+        echo "$scriptname: '$1' doesn't take any values"
+    else
+        echo "$scriptname: '$1' (or '$2') doesn't take any values"
+    fi
+    exit 1
+}
+
+check_v() {
+    if [ x"$2" = x ]; then
+        cat <<EOF
+$scriptname: option '$1' requires an argument. Try '$scriptname --help' for 
more information
+EOF
+        exit 1;
+    fi
+}
+
+
+
+
+# Read the command-line options
+# Separate command-line arguments from options and put the option values
+# into the respective variables.
+#
+# OPTIONS WITH A VALUE:
+#
+#   Each option has three lines because we take into account the three common
+#   formats:
+#   For long option names, '--longname value' and '--longname=value'.
+#   For short option names, '-l value', '-l=value' and '-lvalue'
+#   (where '-l' is the short version of the hypothetical '--longname option').
+#
+#   The first case (with a space between the name and value) is two
+#   command-line arguments. So, we'll need to shift it twice. The
+#   latter two cases are a single command-line argument, so we just need to
+#   "shift" the counter by one.
+#
+#   IMPORTANT NOTE: the ORDER OF THE LATTER TWO cases matters: '-h*' should be
+#   checked only when we are sure that its not '-h=*').
+#
+# OPTIONS WITH NO VALUE (ON-OFF OPTIONS)
+#
+#   For these, we just want the forms of '--longname' or '-l'. Nothing
+#   else. So if an equal sign is given we should definitely crash and also,
+#   if a value is appended to the short format it should crash. So in the
+#   second test for these ('-l*') will account for both the case where we
+#   have an equal sign and where we don't.
+while [ $# -gt 0 ]
+do
+    case "$1" in
+
+    # Input parameters.
+    -h|--hdu)        hdu="$2";                                    check_v "$1" 
"$hdu";    shift;shift;;
+    -h=*|--hdu=*)    hdu="${1#*=}";                               check_v "$1" 
"$hdu";    shift;;
+    -h*)             hdu=$(echo "$1" | sed -e's/-h//');           check_v "$1" 
"$hdu";    shift;;
+    -i|--img)        img="$2";                                    check_v "$1" 
"$img";    shift;shift;;
+    -i=*|--img=*)    img="${1#*=}";                               check_v "$1" 
"$img";    shift;;
+    -i*)             img=$(echo "$1" | sed -e's/-i//');           check_v "$1" 
"$img";    shift;;
+    -H|--imghdu)     imghdu="$2";                                 check_v "$1" 
"$imghdu"; shift;shift;;
+    -H=*|--imghdu=*) imghdu="${1#*=}";                            check_v "$1" 
"$imghdu"; shift;;
+    -H*)             imghdu=$(echo "$1" | sed -e's/-H//');        check_v "$1" 
"$imghdu"; shift;;
+    -r|--racol)      racol="$2";                                  check_v "$1" 
"$racol";  shift;shift;;
+    -r=*|--racol=*)  racol="${1#*=}";                             check_v "$1" 
"$racol";  shift;;
+    -r*)             racol=$(echo "$1" | sed -e's/-r//');         check_v "$1" 
"$racol";  shift;;
+    -d|--deccol)     deccol="$2";                                 check_v "$1" 
"$deccol"; shift;shift;;
+    -d=*|--deccol=*) deccol="${1#*=}";                            check_v "$1" 
"$deccol"; shift;;
+    -d*)             deccol=$(echo "$1" | sed -e's/-d//');        check_v "$1" 
"$deccol"; shift;;
+    -C|--center)     center="$2";                                 check_v "$1" 
"$center"; shift;shift;;
+    -C=*|--center=*) center="${1#*=}";                            check_v "$1" 
"$center"; shift;;
+    -C*)             center=$(echo "$1" | sed -e's/-C//');        check_v "$1" 
"$center"; shift;;
+    -w|--width)      width="$2";                                  check_v "$1" 
"$width";  shift;shift;;
+    -w=*|--width=*)  width="${1#*=}";                             check_v "$1" 
"$width";  shift;;
+    -w*)             width=$(echo "$1" | sed -e's/-w//');         check_v "$1" 
"$width";  shift;;
+    --widthinpix)    widthinpix=1;                                             
           shift;;
+    --widthinpix=*)  on_off_option_error --quiet -q;;
+
+   # Output parameters
+   -k|--keeptmp)           keeptmp=1; shift;;
+   -k*|--keeptmp=*)        on_off_option_error --keeptmp -k;;
+   -t|--tmpdir)            tmpdir="$2";                                  
check_v "$1" "$tmpdir";  shift;shift;;
+   -t=*|--tmpdir=*)        tmpdir="${1#*=}";                             
check_v "$1" "$tmpdir";  shift;;
+   -t*)                    tmpdir=$(echo "$1" | sed -e's/-t//');         
check_v "$1" "$tmpdir";  shift;;
+   -o|--output)            output="$2";                                  
check_v "$1" "$output"; shift;shift;;
+   -o=*|--output=*)        output="${1#*=}";                             
check_v "$1" "$output"; shift;;
+   -o*)                    output=$(echo "$1" | sed -e's/-o//');         
check_v "$1" "$output"; shift;;
+
+    # Non-operating options.
+    -q|--quiet)             quiet=" -q"; shift;;
+    -q*|--quiet=*)          on_off_option_error --quiet -q;;
+    -?|--help)              print_help; exit 0;;
+    -'?'*|--help=*)         on_off_option_error --help -?;;
+    -V|--version)           print_version; exit 0;;
+    -V*|--version=*)        on_off_option_error --version -V;;
+    --cite)                 print_citation; exit 0;;
+    --cite=*)               on_off_option_error --cite;;
+
+    # Unrecognized option:
+    -*) echo "$scriptname: unknown option '$1'"; exit 1;;
+
+    # Not an option (not starting with a `-'): assumed to be the name
+    # of the positions catalog.
+    *)
+      if [ x"$cat" = x ]; then cat="$1";
+      else
+         echo "$scriptname: only one input (position catalog) is expected";
+      fi
+      shift;;
+  esac
+done
+
+
+
+
+
+# Necessary options.
+if [ x"$cat" = x ]; then
+    echo "$scriptname: no argument (position catalog) given"; exit 1
+fi
+if [ x"$img" = x ]; then
+    echo "$scriptname: no refernece image given to '--img'"; exit 1
+fi
+if [ x"$width" = x ]; then
+    echo "$scriptname: no stack width given to '--width'"; exit 1
+fi
+if [ x"$center" = x ]; then
+    echo "$scriptname: no stack center given to '--center'"; exit 1
+fi
+
+
+
+
+# Set the temporary directory name and delete its contents if there is
+# any.
+if [ x"$tmpdir" = x ]; then
+    tmpdir=dither-visualize-tmp
+fi
+if [ -d $tmpdir ]; then rm -f $tmpdir/*
+else                    mkdir $tmpdir
+fi
+
+
+
+
+
+# Simulate the image for each dither pattern.
+counter=1
+elist=$tmpdir/exp-list.txt
+const=$tmpdir/constant.fits
+icopy=$tmpdir/input-copy.fits
+rm -f $elist
+if [ $widthinpix = 0 ]; then wip=""; else wip="--widthinpix"; fi
+asttable $cat -c$racol,$deccol \
+    | while read r d; do
+
+        # Copy the input image HDU into the temporary directory for
+        # modification.
+        astfits $img --copy=$imghdu --output=$icopy
+       astfits $icopy --update=CRVAL1,$r --update=CRVAL2,$d
+
+       # Set all the pixel values to 1.
+       astarithmetic $icopy 1 uint8 constant --output=$const
+       rm $icopy
+
+       # Warp the image to the final pixel grid.
+       warped=$tmpdir/warped-$counter.fits
+       astwarp $const --center=$center --width=$width $wip \
+               --output=$warped
+       rm $const
+
+       # In this scenario, we do not care about the values, so set
+       # all the non-blank pixels to 1.
+       exp=$tmpdir/exp-$counter.fits
+       astarithmetic $warped isnotblank -o$exp
+       rm $warped
+
+       # Add this image to the list of images to add.
+       printf "%s " $exp >> $elist
+
+       # Increment the counter.
+        counter=$((counter+1))
+done
+
+# Stack all the exposure maps.
+astarithmetic $(cat $elist) $(wc $elist -w | awk '{print $1}') \
+             sum -g1 --output=$output
+
+# Delete the temporary directory if necessary
+if [ "$keeptmp" = 0 ]; then
+   rm -rf $tmpdir
+fi