[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 8089ea3 27/32: astscript-radial-profile: Samae
From: |
Mohammad Akhlaghi |
Subject: |
[gnuastro-commits] master 8089ea3 27/32: astscript-radial-profile: Samaeh's changes included, conflicts fixed |
Date: |
Wed, 24 Feb 2021 22:36:20 -0500 (EST) |
branch: master
commit 8089ea3c1bdba27d810c3e5926f3a69caaf56709
Merge: e41c2e4 ad7fb58
Author: Raul Infante-Sainz <infantesainz@gmail.com>
Commit: Raul Infante-Sainz <infantesainz@gmail.com>
astscript-radial-profile: Samaeh's changes included, conflicts fixed
With this commit, I am including the modifications and changes done by
Samaeh into the radial profile script.
---
bin/script/radial-profile.in | 277 +++++++++++++++++++++++++++++++------------
1 file changed, 199 insertions(+), 78 deletions(-)
diff --git a/bin/script/radial-profile.in b/bin/script/radial-profile.in
index b40515c..f1f878d 100644
--- a/bin/script/radial-profile.in
+++ b/bin/script/radial-profile.in
@@ -7,6 +7,7 @@
# Copyright (C) 2020 Raul Infante-Sainz <infantesainz@gmail.com>
# Contributing author(s):
# Copyright (C) 2020 Mohammad Akhlaghi <mohammad@akhlaghi.org>
+# Copyright (C) 2020 Zahra Sharbaf <zahra.sharbaf2@gmail.com>
# Copyright (C) 2020, Free Software Foundation, Inc.
#
# Gnuastro is free software: you can redistribute it and/or modify it under
@@ -32,7 +33,7 @@ set -e
# Default option values (can be changed with options on the command-line).
hdu=1
-rmax=10
+rmax=max
mode=img
x=center
y=center
@@ -51,6 +52,7 @@ j="v"
l="S/N"
quiet=0
prefix=./
+tmpdir=""
output="default"
version=@VERSION@
scriptname=@SCRIPT_NAME@
@@ -253,9 +255,12 @@ do
-l=*|--lname=*) l="${1#*=}"; check_v "$1" "$k";
shift;;
-l*) l=$(echo "$1" | sed -e's/-l//'); check_v "$1" "$k";
shift;;
+
# Output parameters
-k|--keeptemp) k=0; shift;;
-k*|--keeptemp=*) on_off_option_error --keeptemp -k;;
+ --tmpdir) tmpdir="$2"; check_v "$1"
"$tmpdir"; shift;shift;;
+ --tmpdir=*) tmpdir="${1#*=}"; 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;;
@@ -293,34 +298,91 @@ fi
+# If one of X or Y are given the other also needs to be given.
+if [ "z$x" = zcenter ]; then
+ if ! [ "z$y" = zcenter ]; then
+ echo "Center position's Y axis value is given, but not X!"
+ exit 1
+ fi
+else
+ if [ "z$y" = zcenter ]; then
+ echo "Center position's X axis value is given, but not Y!"
+ exit 1
+ fi
+fi
+
+
+
+
+
+# Convert center to image coordinates if necessary
+# ------------------------------------------------
+#
+# If the user gave the central position, and has said its in WCS, then
+# convert them to image mode so we can safely assume image coordianates
+# from now on.
+if ! [ "z$x" = zcenter ]; then
+ if [ $mode = wcs ]; then
+ xy=$(echo "$x $y" \
+ | asttable -c'arith $1 $2 wcstoimg' \
+ --wcsfile=$inputs --wcshdu=$hdu)
+ x=$(echo $xy | awk '{print $1}');
+ y=$(echo $xy | awk '{print $2}');
+ fi
+fi
+
+
+
+
+
-# Calculate the center of the image
-# ---------------------------------
+# Set default central position
+# ----------------------------
#
# If the user don't set the x and y coordinates to be "center" (The
# coordinates of the object), then compute the center of the image for
# constructing the profiles.
-if [ "z$x" = zcenter ] && [ "z$y" = zcenter ]; then
- xpix=$(astfits $inputs --hdu=$hdu | grep NAXIS1 | awk '{print $3}')
- x=$(echo "$(seq $xpix)" | aststatistics --median)
- ypix=$(astfits $inputs --hdu=$hdu | grep NAXIS2 | awk '{print $3}')
- y=$(echo "$(seq $ypix)" | aststatistics --median)
+#
+# In the FITS standard, pixels are counted from 1, and the integers are in
+# the center of the pixel. So after dividing the pixel size of the image by
+# 2, we should add it with 0.5 to be the "center" of the image.
+if [ "z$x" = zcenter ]; then
+ x=$(astfits $inputs --hdu=$hdu | awk '/^NAXIS1/{print $3/2+0.5}')
+ y=$(astfits $inputs --hdu=$hdu | awk '/^NAXIS2/{print $3/2+0.5}')
fi
+
# Calculate the maximum radius
# ----------------------------
#
# If the user set the rmax parameter to "max", then compute the maximum
-# radius possible on the image. Assuming that the object is in one corner of
-# the image, this radius is given by the diagonal of the image
+# radius possible on the image.
+#
+# If the user hasn't given any maximum radius, we give the most reliable
+# maximum radius (where the full circumference will be within the
+# image). If the radius goes outside the image, then the measurements and
+# calculations can be biased, so when the user hasn't provided any maximum
+# radius, we should only confine ourselves to a radius where the results
+# are reliable.
+#
+# Y--------------
+# | | The maximum radius (to ensure the profile
+# y |........* | lies within the image) is the smallest
+# | . | one of these values:
+# | . | x, y, X-x, Y-y
+# --------------
+# 0 x X
+#
if [ "z$rmax" = zmax ]; then
- xpix=$(astfits $inputs --hdu=$hdu | grep NAXIS1 | awk '{print $3}')
- ypix=$(astfits $inputs --hdu=$hdu | grep NAXIS2 | awk '{print $3}')
- rmax=$(astarithmetic "$xpix"f 2f pow "$ypix"f 2f pow + sqrt --quiet)
+ rmax=$(astfits $inputs --hdu=$hdu \
+ | awk '/^NAXIS1/{X=$3} /^NAXIS2/{Y=$3} \
+ END{ x='$x'; y='$y'; \
+ printf("%s\n%s\n%s\n%s", x, y, X-x, Y-y); }' \
+ | aststatistics --minimum )
fi
@@ -334,23 +396,42 @@ fi
# If the user has defined a specific path/name for the output, it will be
# used for saving the output file. If the user do not specify a output name,
# then a default value containing the center and mode will be generated.
-if [ z$output = zdefault ]; then
- bname_ext=$(basename $inputs)
- bname=$(echo $bname_ext | sed
-e"s|.fits|_radial_profile_"$mode"_$x-$y.fits|g")
- output=$prefix$bname
-fi
+bname_prefix=$(basename $inputs | sed 's/\.fits/ /' | awk '{print $1}')
+defaultname=$(pwd)/"$bname_prefix"_radial_profile_$mode"_$x"_"$y"
+if [ z$tmpdir = z ]; then tmpdir=$defaultname
+else tmpdir=$(realpath $tmpdir); fi
+if ! [ -d $tmpdir ]; then mkdir $tmpdir; fi
+if [ z$output = zdefault ]; then output="$defaultname.fits"; fi
-# Generate the aperture parameters
-# --------------------------------
+# Crop image
+# ----------
#
-# Here, an ascii file is generated with the parameters that will be feed
-# into MakeProfiles in order to construct the apertures image.
-aperturestxt=$(echo $output | sed -e"s|.fits|_apertures.txt|g")
-echo "$rmax $x $y 7 $rmax 1 $p $Q $rmax 1" > $aperturestxt
+# Crop the input image around the desired point so we can continue
+# processing only on those pixels (we don't need the other pixels).
+#
+# Crop's output always has the range of pixels from the original image used
+# in the 'ICF1PIX' keyword value. So to find the new center (important if
+# it has sub-pixel positions), we can simply get the first and third value
+# of that string, and convert to the cropped coordinate system. Note that
+# because FITS pixel couting starts from 1, we need to subtract '1'.
+crop=$tmpdir/crop.fits
+cropwidth=$(echo $rmax | awk '{print $1*2+1}')
+astcrop $inputs --hdu=$hdu --center=$x,$y --mode=img \
+ --width=$cropwidth --output=$crop
+dxy=$(astfits $crop -h1 \
+ | grep ICF1PIX \
+ | sed -e"s/'/ /g" -e's/\:/ /g' -e's/,/ /' \
+ | awk '{print $3-1, $5-1}')
+echo "x:$x"
+echo "y:$y"
+echo "cropwidth: $cropwidth"
+echo "dxy: $dxy"
+x=$(echo "$x $cropwidth $dxy" | awk '{if($1>int($2/2)) print $1-$3; else print
int($2/2)+$1-int($1)}')
+y=$(echo "$y $cropwidth $dxy" | awk '{if($1>int($2/2)) print $1-$4; else print
int($2/2)+$1-int($1)}')
@@ -361,24 +442,14 @@ echo "$rmax $x $y 7 $rmax 1 $p $Q $rmax 1" > $aperturestxt
#
# The apertures image is genrated using MakeProfiles with the parameters
# previously specified in the ascii file.
-aperturesfits=$(echo $output | sed -e"s|.fits|_apertures.fits|g")
-astmkprof $aperturestxt --background=$inputs --backhdu=$hdu \
- --mforflatpix --mode=$mode --clearcanvas --type=int16 \
- --circumwidth=$w --replace --output=$aperturesfits \
- --config=$a
-
-
-
-
-
-# Detection signal
-# ----------------
-#
-# To obtain signal to noise ratio radial profile it is necessary to have the
-# detection map. For now, it uses the default option of NoiseChisel.
-detection=$(echo $output | sed -e"s|.fits|_detected.fits|g")
-astnoisechisel $inputs -o$detection
-
+apertures=$tmpdir/apertures.fits
+echo "x:$x"
+echo "y:$y"
+echo "$rmax $x $y 7 $rmax 1 $p $Q $rmax 1" \
+ | astmkprof --background=$crop --backhdu=1 --mforflatpix \
+ --mode=$mode --clearcanvas --type=int16 \
+ --circumwidth=$w --replace --output=$apertures \
+ --config=$a
@@ -386,16 +457,41 @@ astnoisechisel $inputs -o$detection
# Obtain the radial profile
# -------------------------
#
-# The radial profile is obtained using Catalog. In practice, what is done is to
-# obtain a catalogue using the profile image previously generated (the
-# elliptical apertures) and the original input image for computing the values.
-# Also using of SKY_STD extension of NoiseChisel output for computing signal
-# to noise ratio.
-fprofile=$(echo $output | sed -e"s|.fits|_catalogue.fits|g")
-astmkcatalog $aperturesfits -h1 --valuesfile=$inputs --valueshdu=$hdu \
- --instd=$detection --stdhdu=SKY_STD \
- --config=$c -o$fprofile --sigmaclip=$s \
- --ids --$m --sn
+# The radial profile is obtained using Catalog. In practice, what is done is
+# to obtain a catalogue using the segmentation image previously generated
+# (the elliptical apertures) and the original input image for computing the
+# values.
+if [ $b = 1 ]; then
+ fprofile=$tmpdir/catalog-apertures.fits
+ astmkcatalog $apertures -h1 --valuesfile=$crop --valueshdu=1 \
+ --ids --$m --config=$c -o$fprofile
+else
+ # Detection signal
+ #
+ # If the user set the default value of binning variable($b) except 1,
+ # we need to have signal to noise ratio column to obtain radial
+ # profile. To obtain signal to noise ratio for each radial profile; at
+ # first we have to detect signal of noise.
+ detection=$(echo $output | sed -e"s|.fits|_detected.fits|g")
+ astnoisechisel $crop -o$detection
+
+
+
+
+
+ # Obtain the radial profile
+ # -------------------------
+ #
+ # The radial profile is obtained using Catalog. In practice, what is
+ # done is to obtain a catalogue using the segmentation image previously
+ # generated (the elliptical apertures) and the original input image for
+ # computing the values. Also using of SKY_STD extinction of noisechisel
+ # output for computing signal to noise ratio.
+ fprofile=$tmpdir/catalog-apertures.fits
+ astmkcatalog $apertures -h1 --valuesfile=$crop --valueshdu=$hdu \
+ --instd=$detection --stdhdu=SKY_STD --ids --$m --sn \
+ --config=$c -o$fprofile
+fi
@@ -404,14 +500,16 @@ astmkcatalog $aperturesfits -h1 --valuesfile=$inputs
--valueshdu=$hdu \
# Binning data
# ------------
#
-# In order to increase the signal-to-noise ratio of the radial profile, it is
-# possible to bin the data. If the user specify the binning parameter (`-b' or
-# `--binning') different of 1. To do the binning of the data, a small Awk
-# script is used. Because of the behavior of signal to noise ratio is not
-# linear, a different way of binning the S/N column is used. Since the Awk
-# script will print the columns as float values, it is necessary to change the
-# headers. To do that, Sed is used to replace all ocurrences of i32 to f32.
-bprofile=$(echo $output | sed -e"s|.fits|_binned.fits|g")
+# In order to increase the signal-to-noise ratio of the radial profile, it
+# is possible to bin the data. It is done in any case, if the user set the
+# default value of binning variable($b) except 1, because the binning will
+# be equal to 1 so the output binned will be the same as the input. To do
+# the binning of the data, a small Awk script is used.Because of the
+# behavior of signal to noise ratio isn't linear. We used the different
+# way to do binning the S/N column. Since the Awk script will print the
+# columns as float values, it is necessary to change the headers. To do
+# that, Sed is used to replace all ocurrences of i32 to f32.
+bprofile=$tmpdir/catalog-apertures-binned.fits
if [ $b != 1 ]; then
asttable $fprofile \
@@ -476,17 +574,30 @@ if [ z$j = zv ]; then
ycolname=$j$m
fi
-aprofile=$(echo $output | sed -e"s|.fits|_arith.fits|g")
+aprofile=$tmpdir/catalog-arith.txt
echo "# Column 1: $xcolname [modified,f32,] Binned $b and arith $X" > $aprofile
echo "# Column 2: $ycolname [modified,f32,] Binned $b and arith $Y" >>
$aprofile
-echo "# Column 3: $zcolname [modified,f32,] Binned $b and arith $Z" >>
$aprofile
+#
+if [ $b != 1 ]; then
+ echo "# Column 3: $zcolname [modified,f32,] Binned $b and arith $Z" >>
$aprofile
+
+ # Make the appropiate operation over the two columns of the radial profile
+ # and add it to the temporal file with the meta-data.
+ asttable $bprofile \
+ -c'arith $1 '"$X" \
+ -c'arith $2 '"$Y" \
+ -c'arith $3 '"$Z" >> $aprofile
+else
+ # Make the appropiate operation over the two columns of the radial profile
+ # and add it to the temporal file with the meta-data.
+ asttable $bprofile \
+ -c'arith $1 '"$X" \
+ -c'arith $2 '"$Y" >> $aprofile
+fi
+
+
+
-# Make the appropiate operation over the two columns of the radial profile
-# and add it to the temporal file with the meta-data.
-asttable $bprofile \
- -c'arith $1 '"$X" \
- -c'arith $2 '"$Y" \
- -c'arith $3 '"$Z" >> $aprofile
# Finally, read the temporal file with the metadata information as well as
# the radial profile values, and save it as a fits table using Table.
@@ -499,12 +610,20 @@ cat $aprofile | asttable -o$output
# Calculate the Half-Light Radii
# ------------------------------
#
+# This parameter can calulate for that time the user set the default value
+# of binning variable($b) except 1.
+#
# The half-light or 'effective' radius as the radius within which half of
# the object's luminosity is contained.
-halfoflight=$(asttable $output | awk 'NR>=3{if(max<$2){max=$2}}END{print
max/2}')
-radius=$(asttable $output \
- | awk '{if('$(echo $halfoflight)'<=$2){line=$1;light=$2}}END{print
line}')
-echo "Half-Light-Radii: $radius pixels from the center of the object"
+if [ $b != 1 ]; then
+ halfoflight=$(asttable $output | awk 'NR>=3{if(max<$2){max=$2}}END{print
max/2}')
+ radius=$(asttable $output \
+ | awk '{if('$(echo
$halfoflight)'<=$2){line=$1;light=$2}}END{print line}')
+ echo "Half-Light-Radii: $radius pixels from the center of the object"
+fi
+
+
+
@@ -514,11 +633,13 @@ echo "Half-Light-Radii: $radius pixels from the center of
the object"
# ---------------------
#
# If the user has specified this option, temporal files will be removed.
+k=0
if [ $k = 1 ]; then
- rm $aprofile \
- $bprofile \
- $fprofile \
- $detection \
- $aperturestxt \
- $aperturesfits
+ rm $crop \
+ $aprofile \
+ $bprofile \
+ $fprofile \
+ $detection \
+ $apertures
+ rm -df $tmpdir
fi
- [gnuastro-commits] master 4f3d70d 20/32: Binning data in some case, (continued)
- [gnuastro-commits] master 4f3d70d 20/32: Binning data in some case, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master b36c21f 17/32: astscript-radial-profile: modification to use $ instead of c in Table, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master e24caa6 07/32: astscript-radial-profile: user can now specify output column names, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master e1ded16 11/32: Book: minor typos and corrections in radial profile script examples, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 87d68b0 12/32: astscript-radial-profile: removing --rmin option because it is not used, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master aec3057 13/32: Book: describing most important options of radial profile script, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 1b505d1 21/32: Binning the S/N column, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 4294506 32/32: astscript-radial-profile: polished the script and book, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master ad7fb58 26/32: Central coordinates of the cropped image, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master d085379 22/32: Radial-profile script documentation, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 8089ea3 27/32: astscript-radial-profile: Samaeh's changes included, conflicts fixed,
Mohammad Akhlaghi <=
- [gnuastro-commits] master 04a8bec 29/32: Book: simplifying documentation of astscript-radial-profile, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 362819c 30/32: astscript-radial-profile: fixed bug when creating the temporal directory, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master c9647b5 31/32: Imported work on the radial profile script, minor conflict fixed, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 1072b2b 02/32: astscript-radial-profile: added the option -b for binning the data, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 87f5383 23/32: astscript-radial-profile: minor corrections of previous commits, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master e41c2e4 24/32: astscript-radial-profile: sigma clip parameters for measurements added, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 337e46d 18/32: The half-light radii, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 2df5a22 19/32: Signal to noise ratio for each radius, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 846fc39 25/32: Radial profile script modifications, Mohammad Akhlaghi, 2021/02/24
- [gnuastro-commits] master 8e5ab9e 28/32: astscript-radial-profile: simplifying the script, Mohammad Akhlaghi, 2021/02/24