[gnuastro-commits] master 31f0373a 05/69: PSF model: new script for putting a PSF image into a sky position

[Mohammad Akhlaghi]

branch: master
commit 31f0373a8940172b87a7d842873ff9b280322d24
Author: Raul Infante-Sainz <infantesainz@gmail.com>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    PSF model: new script for putting a PSF image into a sky position
    
    With this commit, a new script has been added. It is desired to model the
    scattered light field (from point-like sources) of an astronomical image.
    To do that, it reads an input image, a PSF image, a center coordinates, and
    a flux factor. With these inputs, the script will situate the PSF image
    (scaled by the flux factor) into the center position provided.
    
    Necessary offsets and cropping factors are computed internally. The final
    output is the scaled PSF at the specified position with the same size than
    the original input image. Note that this script works on invididual stars
    and the flux factor has to be provided. As a consequence, to compute an
    entire scattered light field model from several stars, it is necessary to
    sum the different outputs of this script.
---
 bin/script/psf-model-scattered-light.in | 487 ++++++++++++++++++++++++++++++++
 1 file changed, 487 insertions(+)

diff --git a/bin/script/psf-model-scattered-light.in 
b/bin/script/psf-model-scattered-light.in
new file mode 100755
index 00000000..90234612
--- /dev/null
+++ b/bin/script/psf-model-scattered-light.in
@@ -0,0 +1,487 @@
+#!/bin/sh
+
+# Situate a PSF stamp into a given position with the same output image size
+# than the original image. This script is useful to generate a scattered
+# light field image. Run with `--help', or see description under
+# `print_help' (below) for more.
+#
+# Original author:
+#   Raul Infante-Sainz <infantesainz@gmail.com>
+# Contributing author(s):
+#   Mohammad Akhlaghi <mohammad@akhlaghi.org>
+#   Zahra Sharbaf <zahra.sharbaf2@gmail.com>
+# Copyright (C) 2021, Free Software Foundation, Inc.
+#
+# Gnuastro 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.
+#
+# Gnuastro 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 Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+
+
+# Exit the script in the case of failure
+set -e
+
+
+
+
+
+# Default option values (can be changed with options on the command-line).
+hdu=1
+psf=""
+psfhdu=1
+quiet=""
+center=""
+keeptmp=0
+output=""
+tmpdir=""
+mode="img"
+fluxfactor=""
+version=@VERSION@
+scriptname=@SCRIPT_NAME@
+
+
+
+
+
+# Output of `--usage' and `--help':
+print_usage() {
+    cat <<EOF
+$scriptname: run with '--help' for list of options
+EOF
+}
+
+print_help() {
+    cat <<EOF
+Usage: $scriptname [OPTION] FITS-files
+
+This script is part of GNU Astronomy Utilities $version.
+
+This script will consider a center position to situate the PSF image into a
+new empty image. It is intendeed to generate several images that can be
+added to construct the scattered light field.
+
+For more information, please run any of the following commands. In
+particular the first contains a very comprehensive explanation of this
+script's invocation: expected input(s), output(s), and a full description
+of all the options.
+
+     Inputs/Outputs and options:           $ info $scriptname
+     Full Gnuastro manual/book:            $ info gnuastro
+
+If you couldn't find your answer in the manual, you can get direct help from
+experienced Gnuastro users and developers. For more information, please run:
+
+     $ info help-gnuastro
+
+$scriptname options:
+ Input:
+  -h, --hdu=STR           HDU/extension of all input FITS files.
+  -p, --psf=STR           PSF FITS image.
+  -P, --psfhdu=STR        HDU/extension of the PSF image.
+  -m, --mode=STR          Coordinates mode ('wcs' or 'img').
+  -c, --center=FLT,FLT    Center coordinates of the object.
+  -f, --fluxfactor=FLT    Factor by which the PSF is multiplied.
+
+ Output:
+  -o, --output            Output table with the radial profile.
+  -t, --tmpdir            Directory to keep temporary files.
+  -k, --keeptmp           Keep temporal/auxiliar files.
+
+ Operating mode:
+  -?, --help              Print this help list.
+      --cite              BibTeX citation for this program.
+  -q, --quiet             Don't print the list.
+  -V, --version           Print program version.
+
+Mandatory or optional arguments to long options are also mandatory or optional
+for any corresponding short options.
+
+GNU Astronomy Utilities home page: http://www.gnu.org/software/gnuastro/
+
+Report bugs to bug-gnuastro@gnu.org.
+EOF
+}
+
+
+
+
+
+# Output of `--version':
+print_version() {
+    cat <<EOF
+$scriptname (GNU Astronomy Utilities) $version
+Copyright (C) 2020-2021, 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 Raul Infante-Sainz
+EOF
+}
+
+
+
+
+
+# Functions to check option values and complain if necessary.
+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
+        echo "$scriptname: option '$1' requires an argument."
+        echo "Try '$scriptname --help' for more information."
+        exit 1;
+    fi
+}
+
+
+
+
+
+# Separate command-line arguments from options. Then put the option
+# value into the respective variable.
+#
+# OPTIONS WITH A VALUE:
+#
+#   Each option has three lines because we want to all common formats: for
+#   long option names: `--longname value' and `--longname=value'. For short
+#   option names we want `-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 two times. 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 two 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;;
+        -p|--psf)            psf="$2";                                  
check_v "$1" "$psf";  shift;shift;;
+        -p=*|--psf=*)        psf="${1#*=}";                             
check_v "$1" "$psf";  shift;;
+        -p*)                 psf=$(echo "$1"  | sed -e's/-p//');        
check_v "$1" "$psf";  shift;;
+        -P|--psfhdu)         psfhdu="$2";                               
check_v "$1" "$psfhdu";  shift;shift;;
+        -P=*|--psfhdu=*)     psfhdu="${1#*=}";                          
check_v "$1" "$psfhdu";  shift;;
+        -P*)                 psfhdu=$(echo "$1"  | sed -e's/-P//');     
check_v "$1" "$psfhdu";  shift;;
+        -O|--mode)           mode="$2";                                 
check_v "$1" "$mode";  shift;shift;;
+        -O=*|--mode=*)       mode="${1#*=}";                            
check_v "$1" "$mode";  shift;;
+        -O*)                 mode=$(echo "$1"  | sed -e's/-O//');       
check_v "$1" "$mode";  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;;
+        -f|--fluxfactor)     fluxfactor="$2";                           
check_v "$1" "$fluxfactor";  shift;shift;;
+        -f=*|--fluxfactor=*) fluxfactor="${1#*=}";                      
check_v "$1" "$fluxfactor";  shift;;
+        -f*)                 fluxfactor=$(echo "$1"  | sed -e's/-f//'); 
check_v "$1" "$fluxfactor";  shift;;
+
+        # 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="--quiet"; 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)           astfits --cite; 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 input FITS
+        # file name.
+        *) inputs="$1 $inputs"; shift;;
+    esac
+
+done
+
+
+
+
+
+# Basic sanity checks
+# ===================
+
+# If an input image is not given at all.
+if [ x"$inputs" = x ]; then
+    echo "$scriptname: no input FITS image files."
+    echo "Run with '--help' for more information on how to run."
+    exit 1
+fi
+
+# If a PSF image (--psf) is not given at all.
+if [ x"$psf" = x ]; then
+    echo "$scriptname: no PSF image provided."
+    echo "$scriptname: a PSF image '--psf' ('-p') should be provided."
+    exit 1
+fi
+
+# If a flux factor (--fluxfactor) is not given at all.
+if [ x"$fluxfactor" = x ]; then
+    echo "$scriptname: no flux factor provided."
+    echo "$scriptname: value to '--fluxfactor' ('-f') should be provided."
+    exit 1
+fi
+
+# If center coordinates (--center) is not given at all.
+if [ x"$center" = x ]; then
+    echo "$scriptname: no center coordinates provided."
+    echo "$scriptname: values to '--center' ('-c') should be provided."
+    exit 1
+else
+    ncenter=$(echo $center | awk 'BEGIN{FS=","}END{print NF}')
+    if [ x$ncenter != x2 ]; then
+        echo "$scriptname: '--center' (or '-c') only take two values, but 
$ncenter were given."
+        exit 1
+    fi
+fi
+
+# Make sure the value to '--mode' is either 'wcs' or 'img'.
+if [ $mode = "wcs" ] || [ $mode = "img" ]; then
+    junk=1
+else
+    echo "$scriptname: wrong value to --mode (-m) provided."
+    echo "$scriptname: value to '--mode' ('-m') should be 'wcs' or 'img'"
+    exit 1
+fi
+
+
+
+
+
+# Center coordinates and default object label
+# -------------------------------------------
+#
+# Obtain the coordinates of the center.
+xcoord=$(echo "$center" | awk 'BEGIN{FS=","} {print $1}')
+ycoord=$(echo "$center" | awk 'BEGIN{FS=","} {print $2}')
+
+# With the center coordinates, generate a specific label for the object
+# consisting in its coordinates.
+objectid="$xcoord"_"$ycoord"
+
+
+
+
+
+# Define a temporal directory and thefinal output file
+# ----------------------------------------------------
+#
+# Construct the temporary directory. If the user does not specify any
+# directory, then a default one with the base name of the input image will
+# be constructed.  If the user set the directory, then make it. This
+# directory will be deleted at the end of the script if the user does not
+# want to keep it (with the `--keeptmp' option).
+
+# The final output stamp is also defined here if the user does not provide
+# an explicit name. If the user has defined a specific path/name for the
+# output, it will be used for saving the output file. If the user does not
+# specify a output name, then a default value containing the center and
+# mode will be generated.
+bname_prefix=$(basename $inputs | sed 's/\.fits/ /' | awk '{print $1}')
+if [ x$tmpdir = x ]; then \
+  tmpdir=$(pwd)/"$bname_prefix"_psf_stamps
+fi
+
+if [ -d $tmpdir ]; then
+  junk=1
+else
+  mkdir -p $tmpdir
+fi
+
+# Output
+if [ x$output = x ]; then
+  output=$tmpdir/psf-stamp-$objectid.fits
+fi
+
+
+
+
+
+# Transform WCS to IMG center coordinates
+# ---------------------------------------
+#
+# If the original coordinates have been given in WCS or celestial units
+# (RA/DEC), then transform them to IMG (pixel). Here, this is done by using
+# the WCS information from the original input image. If the original
+# coordinates were done in IMG, then just use them.
+if [ $mode = wcs ]; then
+  xycenter=$(echo "$xcoord,$ycoord" \
+                  | asttable  --column='arith $1 $2 wcstoimg' \
+                              --wcsfile=$inputs --wcshdu=$hdu $quiet)
+  xcenter=$(echo "$xycenter" | awk '{print $1}')
+  ycenter=$(echo "$xycenter" | awk '{print $2}')
+else
+  xcenter=$xcoord
+  ycenter=$ycoord
+fi
+
+
+
+
+
+# Scale the PSF in flux
+# ---------------------
+#
+# The PSF absolute pixel values are not relevant since they are used to be
+# normalized somehow. Here, the input PSF is scaled (multiplied) by the
+# specified factor (--fluxfactor) in order to appropiately obtain the PSF
+# brightness.
+psffluxscaled=$tmpdir/psf-flux-scaled-$objectid.fits
+astarithmetic $psf --hdu=$psfhdu $fluxfactor x \
+                   --output=$psffluxscaled $quiet
+
+
+
+
+
+# Get the original inputs size
+# ----------------------------
+#
+# In order to have the final PSF with the same size than the input image,
+# it is necessary to compute the original size of the input image.
+xaxis=$(astfits $inputs --hdu=$hdu \
+                | awk '/^NAXIS1/{print $3}')
+yaxis=$(astfits $inputs --hdu=$hdu \
+               | awk '/^NAXIS2/{print $3}')
+
+
+
+
+
+# Get the center of the PSF
+# -------------------------
+#
+# The center of the PSF is assumed to be at the center of the PSF image.
+# So, the center of the PSF is computed here from the size of the PSF image
+# along the two axis (in pixels).
+xpsfaxis=$(astfits $psf --hdu=$psfhdu \
+                   | awk '/^NAXIS1/{print $3}')
+ypsfaxis=$(astfits $psf --hdu=$psfhdu \
+                  | awk '/^NAXIS2/{print $3}')
+
+# To obtain the center of the PSF (in pixels), just divide by 2 the size of
+# the PSF image.
+xpsfcenter=$(astarithmetic $xpsfaxis 2 / --quiet)
+ypsfcenter=$(astarithmetic $ypsfaxis 2 / --quiet)
+
+
+
+
+
+# Translate the PSF image
+# -----------------------
+#
+# In order to allocate the PSF into the center coordinates provided by the
+# user, it is necessary to compute the appropiate offsets along the X and Y
+# axis. After that, the PSF image is warped using that offsets. Note that
+# in order to account for a 1 pixel offset, it is necessary to subtract the
+# value 1.
+xdiff=$(astarithmetic $xcenter $xpsfcenter - 1.0 - --quiet)
+ydiff=$(astarithmetic $ycenter $ypsfcenter - 1.0 - --quiet)
+
+psftranslated=$tmpdir/psf-translated-$objectid.fits
+astwarp $psffluxscaled --translate=$xdiff,$ydiff \
+                       --output=$psftranslated $quiet
+
+
+
+
+
+# Crop the PSF image
+# ------------------
+#
+# Once the PSF has been situated appropiately into the right position, it
+# is necessary to crop it with a sub-pixel precision as well as with the
+# same size than the original input image. Otherwise the subtraction of the
+# PSF or scattered light field model would not be possible. Here, this
+# cropping is done.
+xrange=$(echo "$xdiff $xaxis $xcenter" \
+              | awk '{if($1<0) \
+                       {i=int($1); \
+                        d=-1*i+1; \
+                        min=d; max=$2+d-1;} \
+                      else {min=1; max=$2} \
+                      i=int($3); {min=min+1; max=max+1}
+                      printf "%d:%d", min, max}')
+
+yrange=$(echo "$ydiff $yaxis $ycenter" \
+              | awk '{if($1<0) \
+                       {i=int($1); \
+                        d=-1*i+1; \
+                        min=d; max=$2+d-1;} \
+                      else {min=1; max=$2} \
+                      i=int($3); {min=min+1; max=max+1}
+                      printf "%d:%d", min, max}')
+
+# Once the necessary crop parameters have been computed, use the option
+# '--section' to get the PSF image with the correct size and the center
+# situated where the user has specified.
+psfcropped=$tmpdir/psf-cropped-$objectid.fits
+astcrop $psftranslated --mode=img --section=$xrange,$yrange \
+                       --output=$psfcropped $quiet
+
+
+
+
+
+# Put nan values as zero values
+# -----------------------------
+#
+# In the vast majority of the cases, having zeros in the outer parts of the
+# astronomical images is not a good idea. It is better to have them as nan
+# values. However, for the scattered light field modeling, it is better to
+# have the pixels where there are no flux as zero values. By having them as
+# zero values, it is possible to make the addition of different modeled
+# stars (using the PSF) or the subtraction without having nan values in the
+# final subtracted image.
+astarithmetic $psfcropped --hdu=1 set-i \
+              i i isblank 0 where --output=$output $quiet
+
+
+
+
+
+# Remove temporary files
+# ----------------------
+#
+# If the user does not specify to keep the temporal files with the option
+# `--keeptmp', then remove the whole directory.
+if [ $keeptmp = 0 ]; then
+    rm -r $tmpdir
+fi