[gnuastro-commits] master 6cc5bade: psf-stamp: keep saturated pixels after applying the SN threshold

[Mohammad Akhlaghi]

branch: master
commit 6cc5badebcfbba7c6b133367d0cafe5cbc922ea5
Author: Sepideh Eskandarlou <sepideh.eskandarlou@gmail.com>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    psf-stamp: keep saturated pixels after applying the SN threshold
    
    Until now, when the '--snthresh' option was applied into the saturated
    stars the output was completely NAN. This was because of the central pixels
    being NaN, so the extracted labels (of the center) became 0 and the whole
    image would become NaN.
    
    With this commit, the saturated area is filled with a non-existant value
    (maximum + 1) and then the threshold is applied. After all other regions
    above the threshold have been removed, the NaN values are returned.
    
    This fixes bug #64635.
---
 NEWS                    |  3 +++
 bin/script/psf-stamp.sh | 30 ++++++++++++++++++++++++------
 2 files changed, 27 insertions(+), 6 deletions(-)

diff --git a/NEWS b/NEWS
index 441585d0..9dc72672 100644
--- a/NEWS
+++ b/NEWS
@@ -220,6 +220,9 @@ See the end of the file for license conditions.
   bug #64610: --copykeys of Fits program finds wrong key when the first
               characters of the keyword names match.
   bug #64615: Fits program's '--pixelareaonwcs' only accurate on equator.
+  bug #64635: astscript-psf-stamp's '--snthresh' option produces fully NaN
+              valued image when a saturated star is in the center. Found
+              and fixed by Sepideh Eskandarlou.
 
 
 
diff --git a/bin/script/psf-stamp.sh b/bin/script/psf-stamp.sh
index 03b7abb9..015cd38b 100644
--- a/bin/script/psf-stamp.sh
+++ b/bin/script/psf-stamp.sh
@@ -673,7 +673,16 @@ EOF
                 --center=$xcenter,$ycenter \
                 --width=$xwidthinpix,$ywidthinpix \
                --output=$cropstd $quiet
-       astarithmetic $cropped_masked -h1 set-v \
+
+        # Fill the NAN pixels with maximum value of the image plus one (so
+        # the fill value is not present in the image).
+        fillval=$(astarithmetic $cropped_masked maxvalue 1 + -q)
+        fill=$tmpdir/cropped-masked-fill-nan-pix-$objectid.fits
+        astarithmetic $cropped_masked set-i i i isblank $fillval \
+                      where --output=$fill
+
+        # Apply the threshold and label the regions above the threshold.
+        astarithmetic $fill -h1 set-v \
                      $cropstd -h1        set-s \
                      v s /               set-sn \
                      v sn $snthresh lt \
@@ -682,14 +691,21 @@ EOF
                       all isnotblank 2 connected-components \
                       --output=$croplab
 
-        # Extract the id of the main object.
+        # Extract the label of the central coordinate.
         id=$(astcrop $croplab -h1 --mode=wcs \
                      --center=$xcoord,$ycoord \
                      --widthinpix --width=1 --oneelemstdout -q)
 
-        # Mask all the pixels that does not belong to the main object.
+        # Mask all the pixels that do not belong to the main object.
         msk=$tmpdir/mask-$objectid.fits
-        astarithmetic $cropsnt $croplab $id ne nan where -g1 --output=$msk
+        astarithmetic $cropsnt $croplab $id ne nan where -g1 \
+                      --output=$msk.fits
+
+        # Set all the originally NaN-valued pixels to NaN again.
+        astarithmetic $msk.fits set-i i i $fillval eq nan where \
+                      --output=$msk
+
+        # Clean up.
        mv $msk $cropped_masked
     fi
 
@@ -850,8 +866,10 @@ if [ x"$normradiusmin" != x   -a   x"$normradiusmax" != x 
]; then
     # Select the values within the requested radial range.
     values=$(asttable $radialprofile $quiet \
                       --range=1,$normradiusmin,$normradiusmax)
-    if ! normvalue=$(echo "$values" | aststatistics --column=2 --$normop \
-                                                    $finalsigmaclip -q 2> 
/dev/null); then
+    if ! normvalue=$(echo "$values" \
+                          | aststatistics --column=2 --$normop \
+                                          $finalsigmaclip -q \
+                          2> /dev/null); then
         normvalue=nan
     fi
 else