[gnuastro-commits] master d92c6dd8 7/9: PSF scripts: adding the --segment option for compuing the flux factor

[Mohammad Akhlaghi]

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

    PSF scripts: adding the --segment option for compuing the flux factor
    
    Until now, the computation of the flux was done without using the
    segmentation image. Indeed, two PSF scripts had deprecated options related
    with the masking of contaminant sources.
    
    With this commit, these options have been changed in order to use the
    option --segment, that takes as an input the Segment's output to mask non
    wanted sources. In the same way, the necessary parts of the tutorial and
    descriptions of these options in the book have been added.
---
 bin/script/psf-create-make-stamp.in |  15 ++---
 bin/script/psf-model-flux-factor.in | 123 ++++++++++++++++++++++--------------
 doc/gnuastro.texi                   |  61 ++++++++----------
 3 files changed, 109 insertions(+), 90 deletions(-)

diff --git a/bin/script/psf-create-make-stamp.in 
b/bin/script/psf-create-make-stamp.in
index 21d95e49..8cfb3f16 100644
--- a/bin/script/psf-create-make-stamp.in
+++ b/bin/script/psf-create-make-stamp.in
@@ -47,13 +47,12 @@ output=""
 tmpdir=""
 segment=""
 axisratio=1
-maskskyval=0
 normradii=""
 sigmaclip=""
 stampwidth=""
 normop="median"
 positionangle=0
-corewidth="5,5"
+corewidth="3,3"
 version=@VERSION@
 scriptname=@SCRIPT_NAME@
 
@@ -98,7 +97,8 @@ $scriptname options:
   -W, --stampwidth=INT    Width of the stamp in pixels.
   -n, --normradii=INT,INT Minimum and maximum radii (in pixels)
                           for computing the normalization value.
-  -m, --segment=STR       Output of Segment (with OBJECTS and CLUMPS).
+  -S, --segment=STR       Output of Segment (with OBJECTS and CLUMPS).
+  -w, --corewidth=INT     Area width of the central object in pixels for 
unmasking.
   -N, --normop=STR        Operator for computing the normalization value
                           (mean, sigclip-mean, etc.).
   -Q, --axisratio=FLT     Axis ratio for ellipse maskprofile (A/B).
@@ -208,12 +208,9 @@ do
         -w|--corewidth)      corewidth="$2";                            
check_v "$1" "$corewidth";  shift;shift;;
         -w=*|--corewidth=*)  corewidth="${1#*=}";                       
check_v "$1" "$corewidth";  shift;;
         -w*)                 corewidth=$(echo "$1"  | sed -e's/-w//');  
check_v "$1" "$corewidth";  shift;;
-        -m|--segment)        segment="$2";                              
check_v "$1" "$segment";  shift;shift;;
-        -m=*|--segment=*)    segment="${1#*=}";                         
check_v "$1" "$segment";  shift;;
-        -m*)                 segment=$(echo "$1"  | sed -e's/-m//');    
check_v "$1" "$segment";  shift;;
-        -v|--maskskyval)     maskskyval="$2";                           
check_v "$1" "$maskskyval";  shift;shift;;
-        -v=*|--maskskyval=*) maskskyval="${1#*=}";                      
check_v "$1" "$maskskyval";  shift;;
-        -v*)                 maskskyval=$(echo "$1"  | sed -e's/-v//'); 
check_v "$1" "$maskskyval";  shift;;
+        -S|--segment)        segment="$2";                              
check_v "$1" "$segment";  shift;shift;;
+        -S=*|--segment=*)    segment="${1#*=}";                         
check_v "$1" "$segment";  shift;;
+        -S*)                 segment=$(echo "$1"  | sed -e's/-S//');    
check_v "$1" "$segment";  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;;
diff --git a/bin/script/psf-model-flux-factor.in 
b/bin/script/psf-model-flux-factor.in
index 40468353..5a05d2c9 100644
--- a/bin/script/psf-model-flux-factor.in
+++ b/bin/script/psf-model-flux-factor.in
@@ -33,21 +33,20 @@ set -e
 # Default option values (can be changed with options on the command-line).
 hdu=1
 psf=""
-mask=""
 rmax=""
 quiet=""
 psfhdu=1
 center=""
 keeptmp=0
-maskhdu=1
 output=""
 tmpdir=""
-corewidth=""
+segment=""
 normradii=""
 sigmaclip=""
 stampwidth=""
 psfprofile=""
 psfprofilehdu=1
+corewidth="3,3"
 normop="median"
 psfprofilecol=""
 version=@VERSION@
@@ -102,9 +101,8 @@ $scriptname options:
   -W, --stampwidth=INT    Width of the stamp in pixels.
   -n, --normradii=INT,INT Minimum and maximum radii (in pixels)
                           for computing the flux factor value.
-  -m, --mask=STR          Segmentation image (sky = 0).
-  -M, --maskhdu=STR       HDU/extension of the segmentation image.
   -w, --corewidth=INT     Area width of the central object in pixels for 
unmasking.
+  -S, --segment=STR       Output of Segment (with OBJECTS and CLUMPS).
   -R, --rmax=FLT          Maximum radius for the radial profile (in pixels).
   -N, --normop=STR        Operator for computing the normalization value
                           (mean, sigclip-mean, etc.).
@@ -228,12 +226,9 @@ do
         -w|--corewidth)      corewidth="$2";                            
check_v "$1" "$corewidth";  shift;shift;;
         -w=*|--corewidth=*)  corewidth="${1#*=}";                       
check_v "$1" "$corewidth";  shift;;
         -w*)                 corewidth=$(echo "$1"  | sed -e's/-w//');  
check_v "$1" "$corewidth";  shift;;
-        -m|--mask)           mask="$2";                                 
check_v "$1" "$mask";  shift;shift;;
-        -m=*|--mask=*)       mask="${1#*=}";                            
check_v "$1" "$mask";  shift;;
-        -m*)                 mask=$(echo "$1"  | sed -e's/-m//');       
check_v "$1" "$mask";  shift;;
-        -M|--maskhdu)        maskhdu="$2";                              
check_v "$1" "$maskhdu";  shift;shift;;
-        -M=*|--maskhdu=*)    maskhdu="${1#*=}";                         
check_v "$1" "$maskhdu";  shift;;
-        -M*)                 maskhdu=$(echo "$1"  | sed -e's/-M//');    
check_v "$1" "$maskhdu";  shift;;
+        -S|--segment)        segment="$2";                              
check_v "$1" "$segment";  shift;shift;;
+        -S=*|--segment=*)    segment="${1#*=}";                         
check_v "$1" "$segment";  shift;;
+        -S*)                 segment=$(echo "$1"  | sed -e's/-S//');    
check_v "$1" "$segment";  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;;
@@ -481,6 +476,46 @@ astcrop $inputs --hdu=$hdu --mode=img \
 
 
 
+# Function to find label of desired region
+# ----------------------------------------
+#
+# Given a certain extension ('CLUMPS' or 'OBJECTS'), find the respective
+# label.
+find_central_label () {
+  # Input arguments
+  hdu=$1
+
+  # Set the parameters.
+  case $hdu in
+      CLUMPS) labname=clump;;
+      OBJECTS) labname=object;;
+      *) cat <<EOF
+$scriptname: ERROR: a bug! Please contact us at bug-gnuastro@gnu.org to fix 
the problem. The first argument to 'find_central_label' function is not 
'CLUMPS' or 'OBJECTS'
+EOF
+      ;;
+  esac
+
+  # Crop the labeled image.
+  cropped_core=$tmpdir/cropped-core-$objectid.fits
+  astcrop $segment --hdu=$hdu --mode=img $quiet --width=$corewidth \
+          --center=$xcenter,$ycenter --output=$cropped_core
+  lab=$(astarithmetic $cropped_core unique --quiet)
+  numlab=$(echo "$lab" | awk '{print NF}')
+  if [ $numlab != 1 ]; then
+    cat <<EOF
+$scriptname: ERROR: there is more than one $labname label in the core region 
around the given coordinate (specified by '--corewidth=$corewidth', containing 
$labname labels: $lab). Therefore it is not possible to unambiguously identify 
a single $labname. Please decrease the box size given to '--corewidth'.
+EOF
+    exit 1
+  fi
+
+  # Write the output in a file.
+  echo $lab > $tmpdir/cropped-core-$labname-$objectid.txt
+}
+
+
+
+
+
 # Crop and unlabel the segmentation image
 # ---------------------------------------
 #
@@ -490,46 +525,42 @@ astcrop $inputs --hdu=$hdu --mode=img \
 # process is as follow:
 #   - Crop the original mask image.
 #   - Crop the original mask image to a central region (core), in order to
-#   compute what is the central object id. This is necessary to unmask this
-#   object.
+#     compute what is the central object id. This is necessary to unmask
+#     this object.
 #   - Compute what is the central object value, using the median value.
 #   - In the original cropped mask, convert all pixels belonging to the
-#   central object to zeros. By doing this, the central object becomes as
-#   sky.
+#     central object to zeros. By doing this, the central object becomes as
+#     sky.
 #   - Mask all non zero pixels in the mask image as nan values.
-if [ x"$mask" != x ]; then
-    cropped_mask=$tmpdir/cropped_mask-$objectid.fits
-    astcrop $mask --hdu=$maskhdu --mode=img \
-            --center=$xcenter,$ycenter \
-            --width=$stampwidth --output=$cropped_mask $quiet
-
-  # If the user does not provide a core width, consider it as the maximum
-  # normalization radius value.
-  if [ x"$corewidth" = x ]; then corewidth="$normradiusmax,$normradiusmax"
-  fi
-
-  # Crop just the center of the target star.
-  cropped_core=$tmpdir/cropped_core-$objectid.fits
-  astcrop $mask --hdu=$maskhdu --mode=img \
+if [ x"$segment" != x ]; then
+
+  # Find the object and clump labels of the target.
+  find_central_label CLUMPS
+  find_central_label OBJECTS
+  clab=$(cat $tmpdir/cropped-core-clump-$objectid.txt)
+  olab=$(cat $tmpdir/cropped-core-object-$objectid.txt)
+
+  # Crop the object and clump image to same size as desired stamp.
+  cropclp=$tmpdir/cropped-clumps-$objectid.fits
+  cropobj=$tmpdir/cropped-objects-$objectid.fits
+  astcrop $segment --hdu=OBJECTS --mode=img \
+          --center=$xcenter,$ycenter \
+          --width=$stampwidth --output=$cropobj $quiet
+  astcrop $segment --hdu=CLUMPS --mode=img \
           --center=$xcenter,$ycenter \
-          --width=$corewidth --output=$cropped_core $quiet
-
-  # Compute the label of the central object from the core cropped image.
-  centrallabel=$(aststatistics $cropped_core --median -q)
-
-  # Unlabel the central object.
-  cropped_unlabel=$tmpdir/cropped_unlabel-$objectid.fits
-  astarithmetic $cropped_mask --hdu=1 \
-                $cropped_mask --hdu=1 \
-                $centrallabel eq 0 where --output=$cropped_unlabel $quiet
-
-  # All objects that are not the central are masked here.
-  cropped_masked=$tmpdir/cropped_masked-$objectid.fits
-  astarithmetic $cropped --hdu=1 \
-                $cropped_unlabel --hdu=1 \
-                0 ne nan where --output=$cropped_masked $quiet
+          --width=$stampwidth --output=$cropclp $quiet
+
+  # Mask all the undesired regions.
+  cropped_masked=$tmpdir/cropped-masked-$objectid.fits
+  astarithmetic $cropped --hdu=1 set-i --output=$cropped_masked \
+                $cropobj --hdu=1 set-o \
+                $cropclp --hdu=1 set-c \
+                                       \
+                c o $olab ne 0 where c $clab eq -1 where 0 gt set-cmask \
+                o o $olab eq 0 where set-omask \
+                i omask cmask or nan where
 else
-    cropped_masked=$cropped
+  cropped_masked=$cropped
 fi
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 5d4c6bfa..7c9eb2e8 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -5991,16 +5991,18 @@ Let's see how we can do that.
 We will use the script @file{astscript-psf-model-flux-factor} for computing 
the factor by which the inner part of the PSF has to be multiplied to have the 
same flux at the junction radius.
 This script is more general and it is used also for estimating the flux 
factors when modeling individual stars once the final PSF is constructed (for 
more on this script, see @ref{Invoking astscript-psf-model-flux-factor}).
 The basic parameters are the normalization radii (i.e., the radius ring at 
which the ratio betwen the radial profiles of the outer and inner part is 
done), and the center of the outer part of the PSF.
+You can also provide a segmentation image with a CLUMPS and OBJECTS extensions 
in order to mask contaminant sources as explained above.
 The center of the inner PSF is asumed to be in the center of the image.
 Let's obtain that flux factor with the following command:
 
 @example
 $ astscript-psf-model-flux-factor outer/stack.fits \
            --psf=inner/stack.fits --center=500,500 \
+           --segment=label/67510-seg.fits \
            --mode=img --normradii=15,16
 @end example
 
-The output of this script is a value (98.8871), and the interpretation of the 
value is the following: it is the number by which you have to multiply the 
inner stack to have the same pixel values (on average) than the outer stack in 
the ring of pixels defined by @option{--normradii=15,16}.
+The output of this script is a value (98.9177), and the interpretation of the 
value is the following: it is the number by which you have to multiply the 
inner stack to have the same pixel values (on average) than the outer stack in 
the ring of pixels defined by @option{--normradii=15,16}.
 Basically, the script has computed the mean of the ratio of the radial 
profiles of the two images at the specified normalization radii ring.
 
 Once we have obtained the flux factor, we are ready to unite the outer and the 
inner part of the PSF.
@@ -6013,7 +6015,7 @@ Since the flux factor was computed for a ring of 1 pixel 
between 15 and 16 pixel
 @example
 $ astscript-psf-create-junction outer/stack.fits \
            --core=inner/stack.fits --radius=15 \
-           --fluxfactor=98.8871 --output=psf-complete.fits
+           --fluxfactor=98.9177 --output=psf-complete.fits
 $ ds9 psf-complete.fits
 @end example
 
@@ -6033,6 +6035,7 @@ $ mkdir junction-radius-test
 $ astscript-psf-model-flux-factor outer/stack.fits \
            --psf=inner/stack.fits --center=250,250 \
            --mode=img --normradii=$RI,$RO \
+           --segment=label/67510-seg.fits \
            --output=junction-radius-test/ffactor-$RI-$RO.txt
 
 $ ffactor=$(cat junction-radius-test/ffactor-$RI-$RO.txt)
@@ -6078,6 +6081,7 @@ $ astscript-psf-model-flux-factor label/67510-seg.fits \
            --mode=wcs \
            --normradii=20,30 \
            --psf=psf-complete.fits \
+           --segment=label/67510-seg.fits \
            --center=202.98395469814,47.234397688023 \
            --output=single-star/fluxfactor.txt
 @end example
@@ -6137,6 +6141,7 @@ $ asttable outer/67510-6-10.fits \
                    --center=$ra,$dec \
                    --normradii=20,30 \
                    --psf=psf-complete.fits \
+                   --segment=label/67510-seg.fits \
                    --output=models/fluxes/flux-$counter.txt
 
               fluxfactor=$(cat models/fluxes/flux-$counter.txt)
@@ -24776,25 +24781,18 @@ This option takes two values separated by a comma 
(@key{,}).
 The first value is the inner radius, the second is the outer radius.
 These two radius define a ring of pixels around the center that is used for 
obtaining the normalization value.
 
-@item -m STR
-@itemx --mask=STR
-Filename of a segmentation image.
-It is possible use a segmentation image to mask all the objects that are not 
the central object.
-The segmentation image specified by this option must have the same dimensions 
than the input image.
-Non object pixels must have pixel values equal to zero, these pixels are not 
masked on the final stamp.
-Object pixels must have pixel values different from zero, these pixels are 
masked with NaN values on the final stamp.
-In order to not mask the central object (in case it is identified in the 
segmentation image as an object), a region around the center is considered from 
the segmentation image.
-The size of that region can be specified by the option @option{--corewidth}, 
see below.
-As a consequence, all pixels with values equal to the central one in the 
segmentation image are not masked in the final stamp image.
-Internally, all of these pixels are converted to zero values so they are not 
masked.
-
-@item -M STR
-@itemx --maskhdu=STR
-The HDU/extension of the segmentation image (@option{--mask}).
+@item -S STR
+@itemx --segment=STR
+Filename of a segmentation image from Segment's output.
+This image is used to mask all objects that are not the central clump/object.
+The script computes an internal mask that combine the objects and clumps to 
mask all objects that do not correspond to the central one.
+All clumps of the central object that are not the central clump are also 
masked.
+The result is that all objects and clumps that contaminate the central source 
are masked.
+The size of the region that is used for identifying the central clump/object 
can be specified by the option @option{--corewidth}, see below.
 
 @item -w INT
 @itemx --corewidth=INT
-Width of the region that is used to compute the central object over the 
segmentation image @option{--mask}, with the goal of not masking it.
+Width of the region that is used to compute the central clump/object over the 
segmentation image @option{--segment}, with the goal of not masking them.
 If a segmentation image is provided, then it is necessary to not mask the 
central object.
 To do that, a central region is considered in order to compute the values of 
the central object on the segmentation image.
 Once it has been identified, that pixels are not masked.
@@ -25055,29 +25053,22 @@ The option takes two values separated by a comma 
(@key{,}).
 The first value is the inner radius, the second is the outer radius.
 These two radius define a ring of pixels around the center that is used for 
obtaining the flux factor value.
 
-@item -m STR
-@itemx --mask=STR
-Filename of the segmentation image.
-It is possible use a segmentation image to mask all the objects that are not 
the central object.
-The segmentation image specified by this option must have the same dimensions 
than the input image.
-Non object pixels must have pixel values equal to zero, these pixels are not 
masked on the final stamp.
-Object pixels must have pixel values different from zero.
-In order to not mask the central object (in case it is identified in the 
segmentation image as an object), a region around the center is considered from 
the segmentation image.
-The size of that region can be specified by the option @option{--corewidth}, 
see below.
-As a consequence, all pixels with values equal to the central one in the 
segmentation image are not masked.
-Internally, all of these pixels are converted to zero values so they are not 
masked.
-
-@item -M STR
-@itemx --maskhdu=STR
-The HDU/extension of the segmentation image (@option{--mask}).
+@item -S STR
+@itemx --segment=STR
+Filename of a segmentation image from Segment's output.
+This image is used to mask all objects that are not the central clump/object.
+The script computes an internal mask that combine the objects and clumps to 
mask all objects that do not correspond to the central one.
+All clumps of the central object that are not the central clump are also 
masked.
+The result is that all objects and clumps that contaminate the central source 
are masked.
+The size of the region that is used for identifying the central clump/object 
can be specified by the option @option{--corewidth}, see below.
 
 @item -w INT
 @itemx --corewidth=INT
-Width of the region that is used to compute the central object over the 
segmentation image @option{--mask}, with the goal of not masking it.
+Width of the region that is used to compute the central clump/object over the 
segmentation image @option{--segment}, with the goal of not masking them.
 If a segmentation image is provided, then it is necessary to not mask the 
central object.
 To do that, a central region is considered in order to compute the values of 
the central object on the segmentation image.
 Once it has been identified, that pixels are not masked.
-With this option, it is possible to control the size of the central region for 
computing the central object values.
+With this option, it is possible to control the size of the central region for 
computing the central object label (on the segmentation image).
 
 @item -R FLT
 @itemx --rmax=FLT