[gnuastro-commits] master 809be22b 67/69: Book: PSF construction tutorial, completed until the outer part

[Mohammad Akhlaghi]

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

    Book: PSF construction tutorial, completed until the outer part
    
    Until now, the usage of 'astscript-psf-create-mask-stamp' script didn't use
    Segment's output, and when a '--mask' was given, it would only use a single
    mask: not benefiting from the fact that Segment provides both clumps and
    objects. Since Segment's clumps and objects are ideal for such scenarios
    and this script is built to run within Gnuastro, and to avoid complicating
    the code, its much easier for the user to simply give the output of Segment
    (when saturated pixels have been properly dealt with; as described in the
    first part of the tutorial).
    
    With this commit, the following changes were made to
    'astscript-psf-create-mask-stamp' (these are not yet documented in the
    respective "Invoking ..." section of the book):
    
     --segment: the new name for '--mask', and it now assumes that the file
       given is the output of Segment.
    
     --maskhdu: removed because we now assume the output of Segment.
    
     --saturation: removed because we now describe how to properly deal with
       saturated pixels in the respective tutorial.
    
     --corewidth: default value set to '5,5'. This is because we now explain
       how to properly deal with saturated pixels, so using a very small region
       will give a more robust estimate of the clump and object IDs for the
       user's given coordinate.
    
     - I noticed that '--mask' was given two times in the output of '--help'!
    
    In relation to these changes, the tutorial was updated until the part where
    the outer PSF is created (using three stars with magnitude range 6 to
    8). Some new sections have also been added that need to be filled before we
    can finish this branch.
---
 bin/script/psf-create-make-stamp.in | 149 +++++++++++++++++++-----------------
 doc/gnuastro.texi                   | 110 ++++++++++++++++++--------
 2 files changed, 158 insertions(+), 101 deletions(-)

diff --git a/bin/script/psf-create-make-stamp.in 
b/bin/script/psf-create-make-stamp.in
index ee3f0212..b1e6892a 100644
--- a/bin/script/psf-create-make-stamp.in
+++ b/bin/script/psf-create-make-stamp.in
@@ -40,22 +40,20 @@ set -e
 
 # Default option values (can be changed with options on the command-line).
 hdu=1
-mask=""
 quiet=""
 center=""
 keeptmp=0
-maskhdu=1
 output=""
 tmpdir=""
+segment=""
 axisratio=1
 maskskyval=0
-corewidth=""
-saturation=0
 normradii=""
 sigmaclip=""
 stampwidth=""
 normop="median"
 positionangle=0
+corewidth="5,5"
 version=@VERSION@
 scriptname=@SCRIPT_NAME@
 
@@ -100,12 +98,7 @@ $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, --mask=STR          Segmentation image (sky = 0).
-  -M, --maskhdu=STR       HDU/extension of the segmentation image.
-  -l, --saturation        saturation limit for masking saturation pixels.
-  -w, --corewidth=INT     Area width of the central object in pixels for 
unmasking.
-  -S, --mask=STR          Segmentation image (sky = 0).
-  -s, --maskhdu=STR       HDU/extension of the segmentation image.
+  -m, --segment=STR       Output of Segment (with OBJECTS and CLUMPS).
   -N, --normop=STR        Operator for computing the normalization value
                           (mean, sigclip-mean, etc.).
   -Q, --axisratio=FLT     Axis ratio for ellipse maskprofile (A/B).
@@ -206,9 +199,6 @@ do
         -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;;
-        -l|--saturation)     saturation="$2";                           
check_v "$1" "$saturation";  shift;shift;;
-        -l=*|--saturation=*) saturation="${1#*=}";                      
check_v "$1" "$saturation";  shift;;
-        -l*)                 saturation=$(echo "$1"  | sed -e's/-l//'); 
check_v "$1" "$saturation";  shift;;
         -n|--normradii)      normradii="$2";                            
check_v "$1" "$normradii";  shift;shift;;
         -n=*|--normradii=*)  normradii="${1#*=}";                       
check_v "$1" "$normradii";  shift;;
         -n*)                 normradii=$(echo "$1"  | sed -e's/-n//');  
check_v "$1" "$normradii";  shift;;
@@ -218,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|--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;;
+        -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;;
@@ -432,6 +419,45 @@ 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
 # ---------------------------------------
 #
@@ -441,58 +467,40 @@ 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
-
-  cropped_core=$tmpdir/cropped_core-$objectid.fits
-  astcrop $mask --hdu=$maskhdu --mode=img \
-                --center=$xcenter,$ycenter \
-                --width=$corewidth --output=$cropped_core $quiet
-
-  # If the cropped core labeled image has sky background values (zero by
-  # default), convert them to nan. This prevent considering the central
-  # object to be the sky if there are a lot of sky background pixels
-  # compared to the labeled region.
-  cropped_core_nonzeros=$tmpdir/cropped_core_nozeros-$objectid.fits
-  astarithmetic $cropped_core set-i $quiet \
-                i i $maskskyval eq nan where --output=$cropped_core_nonzeros
-
-  # Compute the label of the central object from the core cropped image.
-  centrallabel=$(aststatistics $cropped_core_nonzeros --median -q)
-
-  # Check if the user has determined the --saturation option.
-  if [ $saturation = 0 ]; then
-      mask_saturated=""
-  else
-      mask_saturated="i i $saturation gt nan where set-i"
-  fi
-
-  # Unlabel the central object and then mask everything else.
-  cropped_masked=$tmpdir/cropped_masked-$objectid.fits
-  astarithmetic $cropped      --hdu=1 set-i \
-                $cropped_mask --hdu=1 set-m \
-                                            \
-                $mask_saturated  \
-                m m $centrallabel eq $maskskyval where set-a \
-                a $maskskyval ne set-mask  \
-                                 \
-                i mask nan where \
-                --output=$cropped_masked $quiet
+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=$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
 fi
@@ -507,6 +515,7 @@ fi
 # Only if the the user has specified a ring of normalization (--normradii).
 # Otherwise set the normalization value equal to 1.0 (no normalization).
 if [ x"$normradiusmin" != x ] && [ x"$normradiusmax" != x ]; then
+
     # Generate the radial profile of the stamp, since it has been already
     # centered on the center of the object, it is not necessary to give the
     # center coordinates. If the user specifies a maximum radius, use it.
@@ -532,7 +541,8 @@ if [ x"$normradiusmin" != x ] && [ x"$normradiusmax" != x 
]; then
     fi
     normvalue=$(asttable $radialprofile $quiet \
                          --range=1,$normradiusmin,$normradiusmax \
-                         | aststatistics --column=2 --$normop $finalsigmaclip 
-q)
+                         | aststatistics --column=2 --$normop \
+                                         $finalsigmaclip -q)
 else
     normvalue=1.0
 fi
@@ -559,4 +569,3 @@ astarithmetic $cropped_masked --hdu=1 $normvalue / 
--output=$output $quiet
 if [ $keeptmp = 0 ]; then
     rm -r $tmpdir
 fi
-
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 591385b9..fccf084f 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -669,7 +669,10 @@ Tutorial on building the extended PSF
 
 * Preparing input for extended PSF::  Which stars should be used?
 * Saturated pixels and Segment's clumps::  Saturation is a major hurdle!
-* Building outer part of PSF::  Constructing the outer PSF.
+* Building outer part of PSF::  Building the outermost PSF wings.
+* Inner parts of the PSF::      Going towards the PSF center.
+* Uniting the different PSF components::  Merging all the components into one 
PSF.
+* Subtracting the PSF::         Having the PSF, we now want to subtract it.
 
 Library
 
@@ -23498,7 +23501,10 @@ An overview of the process is given in @ref{Overview 
of the PSF scripts}.
 @menu
 * Preparing input for extended PSF::  Which stars should be used?
 * Saturated pixels and Segment's clumps::  Saturation is a major hurdle!
-* Building outer part of PSF::  Constructing the outer PSF.
+* Building outer part of PSF::  Building the outermost PSF wings.
+* Inner parts of the PSF::      Going towards the PSF center.
+* Uniting the different PSF components::  Merging all the components into one 
PSF.
+* Subtracting the PSF::         Having the PSF, we now want to subtract it.
 @end menu
 
 @node Preparing input for extended PSF, Saturated pixels and Segment's clumps, 
Tutorial on building the extended PSF, Tutorial on building the extended PSF
@@ -23635,11 +23641,10 @@ However, there is still another problem: on the edges 
of the vertical ``bleeding
 Such that the pixels just touching the saturated pixels have strong negative 
values.
 These will also cause problems!
 So with the same command, let's dilate the saturated regions by four times.
-We will use 1 and 2 connectivity in sequence to avoid creating a box-like mask.
 
 @example
 $ astarithmetic saturated.fits set-i i i 2500 gt \
-                1 dilate 2 dilate 1 dilate 2 dilate nan where \
+                2 dilate 2 dilate 2 dilate 2 dilate nan where \
                 --output=sat-masked.fits
 $ ds9 sat-masked.fits
 @end example
@@ -23701,8 +23706,8 @@ $ astmkprof --kernel=gaussian,2,5 --oversample=1 \
             -olabel/kernel.fits
 $ astarithmetic flat/image.fits set-i i i 2500 gt \
                 2 dilate 2 dilate 2 dilate 2 dilate nan where \
-                --output=label/sat-masked.fits
-$ astconvolve label/sat-masked.fits --kernel=kernel.fits \
+                --output=flat/sat-masked.fits
+$ astconvolve flat/sat-masked.fits --kernel=kernel.fits \
               --domain=spatial --output=label/sat-masked-conv.fits
 $ astarithmetic label/sat-masked-conv.fits 2 interpolate-maxofregion \
                 --output=label/image-conv.fits
@@ -23718,8 +23723,14 @@ $ ds9 -mecube -zscale label/nc.fits -zoom to fit
 Looking in the @code{SKY} extension of NoiseChisel, we see that there is no 
footprint of the bright stars or of M51, so the detection was good!
 Furthermore, the saturated pixels haven't caused any problems and the central 
clumps of bright stars are now a single clump.
 We can now proceed to estimating the outer PSF in @ref{Building outer part of 
PSF}.
+But first, let's clean up behind our selves (we only need the Segment output)
 
-@node Building outer part of PSF,  , Saturated pixels and Segment's clumps, 
Tutorial on building the extended PSF
+@example
+$ rm label/image-conv.fits label/kernel.fits label/sat-masked.fits \
+     label/sat-masked-conv.fits
+@end example
+
+@node Building outer part of PSF, Inner parts of the PSF, Saturated pixels and 
Segment's clumps, Tutorial on building the extended PSF
 @subsubsection Building outer part of PSF
 In @ref{Preparing input for extended PSF}, we described how to create a 
Segment clump and object map, while accounting for saturated stars.
 So we are now ready to start building the outer parts of the PSF.
@@ -23774,8 +23785,8 @@ Depending on the PSF shape, internal reflections, 
ghosts, saturated pixels, and
 
 The selection of the normalization radii is something that requires a good 
understanding of the data.
 To do that, let's use two useful parameters that will help us in the checking 
of the data: @code{--tmpdir} and @code{--keeptmp}.
-With @code{--tmpdir=checking-normradii} all temporal files, including the 
radial profiles, will be save in that directory.
-With @code{--keeptmp} we won't remove the temporal files, so it is possible to 
have a look at them.
+With @code{--tmpdir=checking-normradii} all temporary files, including the 
radial profiles, will be save in that directory (instead of an 
internally-created name).
+With @code{--keeptmp} we won't remove the temporal files, so it is possible to 
have a look at them (by default the temporary directory gets deleted at the 
end).
 It is necessary to specify the @code{--normradii} even if we don't know yet 
the final values.
 Otherwise the script will not generate the radial profile.
 As a consequence, in this step we put the normalization radii equal to the 
size of the stamps.
@@ -23785,15 +23796,16 @@ In this particular step we set it to 
@code{--normradii=150,160}.
 @example
 $ counter=1
 $ mkdir finding-normradii
-$ asttable stars-mag-8-12.fits \
+$ asttable outer/stars.fits \
            | while read -r ra dec mag; do
-               astscript-psf-create-make-stamp image-masked.fits \
+               astscript-psf-create-make-stamp flat/sat-masked.fits \
                     --mode=wcs \
+                    --stampwidth=500 \
                     --center=$ra,$dec \
-                    --stampwidth=150 \
-                    --normradii=150,160 \
+                    --normradii=250,260 \
+                    --segment=label/seg.fits \
                     --output=finding-normradii/$counter.fits \
-                    --tmpdir=checking-normradii --keeptmp; \
+                    --tmpdir=finding-normradii --keeptmp; \
                counter=$((counter+1)); \
              done
 @end example
@@ -23803,12 +23815,12 @@ For example using Topcat or any other plotting tool.
 In the same way, it is always good to check the different stamps to ensure the 
quality and posible two dimensional features that are difficult to detect from 
the radial profiles (i.e., ghosts, internal reflections, etc.).
 
 @example
-$ topcat checking-normradii/rprofile*.fits
-$ ds9 checking-normradii/cropped*.fits
+$ topcat finding-normradii/rprofile*.fits
+$ ds9 finding-normradii/cropped-masked*.fits
 @end example
 
-After some study of this data, we could say that a good normalization ring is 
those pixels between R=20 and R=30 pixels.
-Such ring ensures having a high number of pixels so the estimation of the flux 
normalization will be robust.
+After some study of this data, we could say that a good normalization ring is 
those pixels between R=50 and R=60 pixels.
+Such a ring ensures having a high number of pixels so the estimation of the 
flux normalization will be robust.
 Also, at such distance from the center the signal to noise is high and there 
are not obvious features that can affect the normalization.
 
 We later need the normalization radii in the next steps also.
@@ -23818,30 +23830,66 @@ Furthermore, since there are several stars (as you 
saw from the output of the pr
 @example
 $ counter=1
 $ IMAGE=image-masked.fits
-$ NORMRADII_INNER=20
-$ NORMRADII_OUTER=30
-$ mkdir stamps-outer
-$ asttable stars-mag-8-12.fits \
+$ NORMRADII_INNER=50
+$ NORMRADII_OUTER=60
+$ mkdir outer/stamps
+$ asttable outer/stars.fits \
            | while read -r ra dec mag; do
-               astscript-psf-create-make-stamp $IMAGE \
+               astscript-psf-create-make-stamp flat/sat-masked.fits \
                     --mode=wcs \
+                    --stampwidth=500 \
                     --center=$ra,$dec \
-                    --stampwidth=150 \
-                    --normradii=$NORMRADII_INNER,$NORMRADII_OUTER \
-                    --output=stamps-outer/$counter.fits; \
+                    --segment=label/seg.fits \
+                    --output=outer/stamps/$counter.fits \
+                    --normradii=$NORMRADII_INNER,$NORMRADII_OUTER; \
                counter=$((counter+1)); \
              done
 @end example
 
-After the stamps are created, they are stacked together with a simple 
Arithmetic command.
-The stack is done using the median operator and the stacked image is the PSF 
we are looking for: @file{psf.fits}.
+After the stamps are created, we need to stack them together with a simple 
Arithmetic command (see @ref{Stacking operators}).
+The stack is done using the median operator and the stacked image is the outer 
PSF that we are looking for: @file{psf-outer.fits}.
+
+@example
+$ astarithmetic outer/stamps/*.fits -g1 \
+          $(ls outer/stamps/*.fits | wc -l) median \
+          --output=outer/stack.fits
+@end example
+
+@noindent
+With the command below, we'll compare this stacked PSF with the images of the 
individual stars that were used to create it.
+You clearly see that the number of masked pixels is significantly decreased 
and the PSF is much more ``cleaner'').
+
+@example
+$ ds9 outer/stack.fits outer/stamps/*.fits
+@end example
+
+However, the bleeding, vertical saturation in the center still remains.
+Also, because we didn't have too many images (only three!), some small regions 
still remain that were (by chance!) masked in all three.
+If we had more bright stars in our selected magnitude range, we could have 
filled those outer remaining patches.
+
+@c -------------------------------------------
+@c ------- This part may be removed ----------
+To fill those regions, we can use the radial profile of the stacked PSF.
+Obtaining the radial profile for a simple (square, with a single, circular 
object), image is very easy:
 
 @example
-$ astarithmetic starstamp*.fits -g1 \
-          $(ls starstamp*.fits | wc -l) median \
-          --output=psf.fits
+$ astscript-radial-profile outer/stack.fits -oouter/stack-radial.fits
+$ echo "1 251 251 8 300 0 0 1 1 1" \
+       | astmkprof --background=outer/stack.fits --clearcanvas \
+                   --customtable=outer/stack-radial.fits
+
 @end example
+@c -------------------------------------------
+
+@node Inner parts of the PSF, Uniting the different PSF components, Building 
outer part of PSF, Tutorial on building the extended PSF
+@subsubsection Inner parts of the PSF
+
+@node Uniting the different PSF components, Subtracting the PSF, Inner parts 
of the PSF, Tutorial on building the extended PSF
+@subsubsection Uniting the different PSF components
+
 
+@node Subtracting the PSF,  , Uniting the different PSF components, Tutorial 
on building the extended PSF
+@subsubsection Subtracting the PSF
 With the commands above PSF has been constructed and now it is going to be 
used for modeling and subtracting each individual star.
 By construction, the pixel values of the PSF came from the normalization of 
the individual stamps.
 As a consequence, when modeling each individual star, it is necessary to 
compute what is the factor by which the PSF has to be scaled.