[gnuastro-commits] master 844fe9f2 08/11: Book: adding new information to the PSF subtraction tutorial

[Mohammad Akhlaghi]

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

    Book: adding new information to the PSF subtraction tutorial
    
    Until now, there was some missing information in the PSF subtraction
    tutorial about the new features added to the PSF scripts.
    
    With this commit, the tutorial describing how to model and subtract the PSF
    has been updated. Now it includes information on how to use the new
    features. Also, a bug in 'psf-scale-factor' that I found when writing the
    tutorial has been corrected. It was about the injection of the fake
    astrometry into the PSF image.
---
 bin/script/psf-scale-factor.sh | 170 +++++++++++++++++----------------
 doc/gnuastro.texi              | 209 +++++++++++++++++++++++++++++++++++++----
 2 files changed, 279 insertions(+), 100 deletions(-)

diff --git a/bin/script/psf-scale-factor.sh b/bin/script/psf-scale-factor.sh
index 62b7b5c5..37c23891 100644
--- a/bin/script/psf-scale-factor.sh
+++ b/bin/script/psf-scale-factor.sh
@@ -535,18 +535,78 @@ ywidth=$(astarithmetic $normradiusmax float32 2.0 x 1.0 + 
--quiet)
 
 
 
-# Warp the PSF image
-# ------------------
+# Transform WCS to IMG center coordinates
+# ---------------------------------------
 #
-# In principle, the PSF does not have WCS information. Here, fake WCS info
-# is included in order to use warp and center/re-size the PSF to the same
-# size than the warped image. Fake WCS is constructed as follow: reference
-# pixel is the center pixel. This pixel is assigned to be RA, Dec = 180.0,
-# 0.0 deg. Pixel scale is set to 1.0 arcsec/pixel in both axis.
-crval1=180.0
-crval2=0.0
-cdelt1=1.0
-cdelt2=1.0
+# 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
+
+    # For the IMAGE: center coordinates in wcs and img
+    xycenter=$(echo "$xcoord,$ycoord" \
+                   | asttable  --column='arith $1 $2 wcs-to-img' \
+                               --wcsfile=$inputs --wcshdu=$hdu $quiet)
+
+    xcwcs=$xcoord
+    ycwcs=$ycoord
+    xcimg=$(echo "$xycenter" | awk '{print $1}')
+    ycimg=$(echo "$xycenter" | awk '{print $2}')
+    inputs_wcs=$inputs
+
+    # For the PSF: center coordinates in wcs and img
+    crval1=$xcwcs
+    crval2=$ycwcs
+    cdelt1=$(astfits $inputs --hdu=$hdu --pixelscale --quiet | awk '{print 
$1}')
+    cdelt2=$(astfits $inputs --hdu=$hdu --pixelscale --quiet | awk '{print 
$2}')
+
+else
+    # In principle, the PSF does not have WCS information. Here, fake WCS info
+    # is included in order to use warp and center/re-size the PSF to the same
+    # size than the warped image. Fake WCS is constructed as follow: reference
+    # pixel is the center pixel. This pixel is assigned to be RA, Dec = 180.0,
+    # 0.0 deg. Pixel scale is set to 1.0 arcsec/pixel in both axis.
+    crval1=180.0
+    crval2=0.0
+    cdelt1=1.0
+    cdelt2=1.0
+    # For the IMAGE: center coordinates in wcs and img
+    naxis1=$(astfits $inputs --hdu=$hdu --keyval NAXIS1 --quiet)
+    naxis2=$(astfits $inputs --hdu=$hdu --keyval NAXIS2 --quiet)
+    center1=$(echo $naxis1 | awk '{print $1/2}')
+    center2=$(echo $naxis2 | awk '{print $1/2}')
+
+    fake_wcs=$tmpdir/fake-wcs-$objectid.fits
+    echo "1 $center1 $center2 1 1 1 1 1 1 1" \
+        | astmkprof --oversample=1 \
+                    --type=uint8 \
+                    --output=$fake_wcs \
+                    --cdelt=$cdelt1,$cdelt2 \
+                    --crval=$crval1,$crval2 \
+                    --crpix=$center1,$center2 \
+                    --cunit="arcsec,arcsec" \
+                    --mergedsize=$naxis1,$naxis2
+
+    # Inject the fake WCS into the original image, its HDU becomes 1.
+    hdu=1
+    inputs_wcs=$tmpdir/input-wcsed-$objectid.fits
+    astarithmetic $inputs --hdu $hdu \
+                  --wcsfile=$fake_wcs --output $inputs_wcs
+
+    # Center coordinates in IMG and WCS
+    xycenter=$(echo "$xcoord,$ycoord" \
+                   | asttable  --column='arith $1 $2 img-to-wcs' \
+                               --wcsfile=$inputs_wcs --wcshdu=1 $quiet)
+
+    xcimg=$xcoord
+    ycimg=$ycoord
+    xcwcs=$(echo "$xycenter" | awk '{print $1}')
+    ycwcs=$(echo "$xycenter" | awk '{print $2}')
+
+fi
+
+# For the PSF
 psfnaxis1=$(astfits $psf --hdu=$psfhdu --keyval NAXIS1 --quiet)
 psfnaxis2=$(astfits $psf --hdu=$psfhdu --keyval NAXIS2 --quiet)
 psfcenter1=$(echo $psfnaxis1 | awk '{print $1/2+0.5}')
@@ -558,7 +618,6 @@ echo "1 $psfcenter1 $psfcenter2 1 1 1 1 1 1 1" \
     | astmkprof --oversample=1 \
                 --type=uint8 \
                 --output=$psf_fake_wcs \
-                --cunit="arcsec,arcsec" \
                 --cdelt=$cdelt1,$cdelt2 \
                 --crval=$crval1,$crval2 \
                 --crpix=$psfcenter1,$psfcenter2 \
@@ -575,14 +634,17 @@ psfxywcs=$(echo "$psfcenter1,$psfcenter2" \
                             --wcsfile=$psf_wcs --wcshdu=1 $quiet \
                             | awk '{printf "%.10e,%.10e\n", $1, $2}')
 
-# Warp the PSF image around its center with the same size than the image
+# Warp the PSF
+# ------------
+#
+# The PSF image is warped around its center
 psf_warped=$tmpdir/psf-warped.fits
 astwarp $psf_wcs \
         --hdu=1 \
         --widthinpix \
         --center=$psfxywcs \
         --width=$xwidth,$ywidth \
-        --output=$psf_warped $quiet \
+        --output=$psf_warped $quiet
 
 
 
@@ -607,64 +669,6 @@ echo "1 $xradcenter $yradcenter 7 $maxradius 0 0 1 1 1" \
 
 
 
-# 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 wcs-to-img' \
-                               --wcsfile=$inputs --wcshdu=$hdu $quiet)
-
-    # Center coordinates in WCS and IMG
-    xcwcs=$xcoord
-    ycwcs=$ycoord
-    xcimg=$(echo "$xycenter" | awk '{print $1}')
-    ycimg=$(echo "$xycenter" | awk '{print $2}')
-    inputs_wcs=$inputs
-
-else
-    # Fake WCS for the image is taken from the already constructed fake WCS
-    # of the PSF image. The only difference is that the crpix correspond to
-    # the center of the image (instead of the PSF image).
-    naxis1=$(astfits $inputs --hdu=$hdu --keyval NAXIS1 --quiet)
-    naxis2=$(astfits $inputs --hdu=$hdu --keyval NAXIS2 --quiet)
-    center1=$(echo $naxis1 | awk '{print $1/2}')
-    center2=$(echo $naxis2 | awk '{print $1/2}')
-
-    fake_wcs=$tmpdir/fake-wcs-$objectid.fits
-    echo "1 $center1 $center2 1 1 1 1 1 1 1" \
-        | astmkprof --oversample=1 \
-                    --type=uint8 \
-                    --output=$fake_wcs \
-                    --cdelt=$cdelt1,$cdelt2 \
-                    --crval=$crval1,$crval2 \
-                    --crpix=$center1,$center2 \
-                    --cunit="arcsec,arcsec" \
-                    --mergedsize=$naxis1,$naxis2
-
-    # Inject the fake WCS into the original image, its HDU becomes 1.
-    hdu=1
-    inputs_wcs=$tmpdir/input-wcsed-$objectid.fits
-    astarithmetic $inputs --hdu $hdu \
-                  --wcsfile=$fake_wcs --output $inputs_wcs
-
-    # Center coordinates in IMG and WCS
-    xycenter=$(echo "$xcoord,$ycoord" \
-                   | asttable  --column='arith $1 $2 img-to-wcs' \
-                               --wcsfile=$inputs_wcs --wcshdu=1 $quiet)
-
-    xcimg=$xcoord
-    ycimg=$ycoord
-    xcwcs=$(echo "$xycenter" | awk '{print $1}')
-    ycwcs=$(echo "$xycenter" | awk '{print $2}')
-
-fi
-
-
 
 
 
@@ -802,25 +806,29 @@ EOF
                                   --wcsfile=$fake_wcs --output=$objects_wcs
                 fi
 
-                # Warp the object and clump images.
-                warped_clp=$tmpdir/warped-clumps-$objectid-$label_shift.fits
-                warped_obj=$tmpdir/warped-objects-$objectid-$label_shift.fits
-                astwarp $clumps_wcs \
+               # Crop the object and clump images. We don't use Warp
+               # because it is not able to properly account integer values
+               # images because they become float32 and thus the
+               # segmentation clumps/objects are not correct.
+                croped_clp=$tmpdir/croped-clumps-$objectid-$label_shift.fits
+                croped_obj=$tmpdir/croped-objects-$objectid-$label_shift.fits
+                astcrop $clumps_wcs \
+                        --mode=wcs \
                         --widthinpix \
                         --width=$xwidth,$ywidth \
-                        --output=$warped_clp $quiet \
+                        --output=$croped_clp $quiet \
                         --center=$xcwcs_shift,$ycwcs_shift
-                astwarp $objects_wcs \
+                astcrop $objects_wcs \
                         --widthinpix \
                         --width=$xwidth,$ywidth \
-                        --output=$warped_obj $quiet \
+                        --output=$croped_obj $quiet \
                         --center=$xcwcs_shift,$ycwcs_shift
 
                 # Mask all the undesired regions.
                 warped_masked=$tmpdir/warped-masked-$objectid-$label_shift.fits
                 astarithmetic $warped     --hdu=1 set-i  \
-                              $warped_obj --hdu=1 set-o \
-                              $warped_clp --hdu=1 set-c \
+                              $croped_obj --hdu=1 set-o \
+                              $croped_clp --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 \
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index ea589c2b..f2c179df 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -6978,23 +6978,31 @@ The sigma-clipped median of these pixels is finally 
reported.
 Unlike the radial profile method (which averages over a circular/elliptical 
annulus for each radius), this method imposes no a-priori shape on the PSF.
 This makes it very useful for complex PSFs (like the case here).
 
-To continue, let's remove the temporary directory and re-run the script but 
with @option{--quiet} mode so we can put the output in a shell variable.
+To continue, let's compute the scaling and keep it into a plain text file.
+Also, remove the temporary directory used previously.
 
 @example
+$ mkdir uniting
 $ rm -r stack_psfmodelscalefactor
-$ scale=$(astscript-psf-scale-factor outer/stack.fits \
-                   --psf=inner/stack.fits --center=501,501 \
-                   --mode=img --normradii=10,15 --quiet)
-$ echo $scale
+$ astscript-psf-scale-factor outer/stack.fits \
+           --psf=inner/stack.fits --center=501,501 \
+           --mode=img --normradii=10,15 --output=uniting/scale.txt --quiet
+$ cat uniting/scale.txt
 @end example
 
+As you can see, into the file there is a one-line header and then three values.
+The first value is the scaling factor that we are looking for.
+The other two values are the coordinates of the center position.
+At this point, we are only interested in the first value (the scaling flux 
factor), the other two correspond to the coordinates of the center and we will 
see them later in the subtraction of the PSF step.
 Now that we know the scaling factor, we are ready to unite the outer and the 
inner part of the PSF.
 To do that, we will use the script @file{astscript-psf-unite} with the command 
below (for more on this script, see @ref{Invoking astscript-psf-unite}).
 The basic parameters are the inner part of the PSF (given to 
@option{--inner}), the inner part's scale factor (@option{--scale}), and the 
junction radius (@option{--radius}).
 The inner part is first scaled, and all the pixels of the outer image within 
the given radius are replaced with the pixels of the inner image.
 Since the flux factor was computed for a ring of pixels between 10 and 15 
pixels, let's set the junction radius to be 12 pixels (roughly in between 10 
and 15):
+Note that we need to avoid the header line (that starts with "#") and get only 
the first column (the scaling factor), to do that we use AWK.
 
 @example
+scale=$(cat uniting/scale.txt | awk '!/^#/{print $1}')
 $ astscript-psf-unite outer/stack.fits \
            --inner=inner/stack.fits --radius=12 \
            --scale=$scale --output=psf.fits
@@ -7018,13 +7026,14 @@ So let's choose a bad normalization radial range (50 to 
60 pixels) and unite the
 The last command will open the two PSFs together in DS9, you should be able to 
immediately see the discontinuity in the union radius.
 
 @example
-$ scale=$(astscript-psf-scale-factor outer/stack.fits \
+$ astscript-psf-scale-factor outer/stack.fits \
                    --psf=inner/stack.fits --center=501,501 \
-                   --mode=img --normradii=50,60 --quiet)
+                   --mode=img --normradii=50,60 --output=uniting/scale-bad.txt 
--quiet
 
+scale_bad=$(cat uniting/scale-bad.txt | awk '!/^#/{print $1}')
 $ astscript-psf-unite outer/stack.fits \
            --inner=inner/stack.fits --radius=55 \
-           --scale=$scale --output=psf-bad.fits
+           --scale=$scale_bad --output=psf-bad.fits
 
 $ astscript-fits-view psf-bad.fits psf.fits --ds9scale=minmax \
            --ds9extra="-scale limits 0 50 -match scale" \
@@ -7073,27 +7082,41 @@ Finally, once the stars have been modeled by the PSF, 
we will subtract it.
 
 First, let's start with a single star.
 Later, when the basic idea has been explained, we will generalize the method 
for any number of stars.
-With the following command we obtain the coordinates (RA and DEC) and 
magnitude of the brightest star in the image (which is on the top edge of the 
image):
+With the following command we obtain the coordinates (RA and DEC) of the 
fainter star that we have used to construct the outer part of the PSF (note the 
@option{--tail 1}).
 
 @example
 $ mkdir single-star
-$ center=$(asttable flat/67510-bright.fits --sort phot_g_mean_mag \
-                    --column=ra,dec --head 1 \
-                    | awk '@{printf "%s,%s", $1, $2@}')
+$ center=$(asttable outer/67510-6-10.fits --sort phot_g_mean_mag \
+                    --column=ra,dec --tail 1 \
+                    | awk '{printf "%.8f,%.8f", $1, $2}')
 $ echo $center
 @end example
 
-With the center position of that star, let's obtain the flux factor using the 
same normalization ring we used for the creation of the outer part of the PSF:
+With the center position of that star, let's obtain the flux factor using the 
same normalization ring we used for the creation of the outer part of the PSF, 
this parameter is saved into a plain text file:
 
 @example
-$ scale=$(astscript-psf-scale-factor label/67510-seg.fits \
-                   --mode=wcs --quiet \
-                   --psf=psf.fits \
-                   --center=$center \
-                   --normradii=10,15 \
-                   --segment=label/67510-seg.fits)
+$ astscript-psf-scale-factor label/67510-seg.fits \
+            --mode=wcs --quiet \
+            --psf=psf.fits \
+            --center=$center \
+            --normradii=10,15 \
+            --segment=label/67510-seg.fits \
+            --output=single-star/scale.txt 
+@end example
+
+Let's have a look at the file generated, it contains three columns.
+The first column is the scaling flux factor.
+The second and third columns correspond to the coordinates of the center 
position.
+Since this file has a first-line header starting with #, we avoid that line 
and read the first column by using AWK.
+
+@example
+$ cat single-star/scale.txt
+$ scale=$(cat single-star/scale.txt | awk '!/^#/@{print $1@}')
 @end example
 
+For now, just ignore the other two columns that correspond to the center 
position since they are the same than the center position we provided with the 
option @option{--center}.
+We will look into them later.
+
 Now we have all the information necessary to model the star using the PSF: the 
position on the sky and the flux factor.
 Let's use this data with the script @file{astscript-psf-subtract} for modeling 
this star and have a look with DS9.
 
@@ -7109,6 +7132,154 @@ $ astscript-fits-view label/67510-seg.fits 
single-star/subtracted.fits \
            --ds9center=$center --ds9mode=wcs --ds9extra="-zoom 4"
 @end example
 
+After opening the two images at the same time, you will realize that the star 
at the center has been modeled and subtracted.
+There are some residuals that remains as par of the natural subtraction 
because the PSF is an averaged model of several stars, while each star has a 
particular shape.
+But so far, with just a few stars we have been able to model this star and 
subtract it quite accurately.
+
+Now, let's see a couple of interesting parameters that we skip in the steps 
above.
+Imagine, for example, that for some reason, the coordinate of this star is not 
accurate (i.e., the center position does not correspond to the center of the 
star on the image).
+There are many reasons why this could happen.
+For example, the center was computed by Gaia (or any other catalog) and the 
image was taken some time after, and because the proper motion of the star, the 
coordinates from the catalog and the coordinate on the image are not the same.
+Other possibility is that the astrometric solution of the image is not 
accurate enough and it has offests.
+To take into account these and other possibilities, the 
@file{astscript-scale-factor} script is able to consider different shifts on 
the center position.
+Let's see this in practice.
+
+First, we are going to modify the coordinates of the star in order to simulate 
a small offset between the center we provide and the real center position of 
the star on the image.
+
+@example
+$ original_center=$center
+$ modified_center=$(echo $center | awk -F',' '@{printf "%.8f,%.8f", 
$1+0.55/3600.0, $2@}')
+
+$ echo $original_center
+$ echo $modified_center
+@end example
+
+Once we have modified by 1 pixel the RA position, let's see how this changes 
the subtraction of the star.
+In this step, we should re-calculate the scaling factor, but since we just 
want to check the effect of the 1-pixel offset, we consider the same flux 
scaling factor already computed.
+
+@example
+$ astscript-psf-subtract label/67510-seg.fits \
+           --mode=wcs \
+           --psf=psf.fits \
+           --scale=$scale \
+           --center=$modified_center \
+           --output=single-star/subtracted-modified-center.fits
+
+$ astscript-fits-view single-star/subtracted.fits 
single-star/subtracted-modified-center.fits \
+           --ds9center=$center --ds9mode=wcs --ds9extra="-zoom 4"
+@end example
+
+With such small shift in the RA coordinate, it is difficult to see any 
difference by eye.
+To help checking the differences, blink between the two frames of the DS9 
image (Frame/Blink).
+Do you see now the changes?
+Which one is a better subtraction?
+It is difficult to say because they are very similar.
+In fact, since this is a real observation, there is not such thing as the 
correct center of the star, but different centers computed using different 
techniques (Gaia, NoiseChisel+Segment+Catalog), etc.
+To select the best value we could shift the center position along the x and y 
directions several times, make the subtraction of all of them, and finally 
check which one gives the "better subtraction".
+This is actually the goal of the options @options{--xshift} and 
@option{--yshifts}.
+
+These two options allow you to consider shifts along the two axis for the 
center position in the computation of the scaling flux factor.
+Let's see this in practice.
+With the command below, the flux scaling factor will be estimated considering 
the original center position and several shifts along the x and y directions.
+The units of the shifts are in pixels, in such a way that 
@option{--xshifts=-1.0,1.0,+1.0} will generate the following shifts along the 
x-axis: -1.0, 0.0, 1.0 pixels.
+With the command below we also ask for the same along the y-axis, so a total 
of 9 shifts will be computed.
+We also request for keepint all temporary files in the @file{checks} directory.
+
+@example
+$ astscript-psf-scale-factor label/67510-seg.fits \
+            --mode=wcs --quiet \
+            --psf=psf.fits \
+            --center=$center \
+            --normradii=10,15 \
+            --xshifts=-1.0,1.0,+1.0 \
+            --yshifts=-1.0,1.0,+1.0 \
+            --segment=label/67510-seg.fits \
+            --output=single-star/scale-shifts.txt \
+            --tmpdir=checks --keeptmp
+@end example
+
+Let's see what happened.
+For each shifted position, a scaling flux factor was computed within the ring 
defined by @option{--normradii}, the standard deviation of such ring was also 
kept.
+Once all of them were computed, the shifted position that gave the smallest 
standard deviation value was selected to be the best one.
+As the final result, in the @file{single-star/scale.txt} there are three 
columns: the scaling factor, and the shifted coordinates.
+It is up to you to use all of only some of them for the final subtraction of 
the PSF.
+
+As an advice, it is usefult to check the interemediate results to get some 
clues in case you obtain some unexpected results.
+This can be done by having a look at the files generated in the @file{checks} 
directory.
+In particular, all ring shifted images can be visually inspectioned.
+The complete results are kept in a @file{fits} table with the following 
columns: 1. best shifted ra coordinate, 2. best shifted dec coordinate, 3. flux 
scaling factor, 4. flux scaling factor standard deviation, 5. ra shift in wcs 
units, 6. dec shift in wcs units, 7. x shift in pixels, 8 y shift in pixels.
+
+@example
+$ astscript-fits-view checks/for-factor-*-xs_ys_*.fits
+
+$ asttable checks/stats-all-*.fits  -Y
+@end example
+
+The final (best) result is saved into the output file 
@file{single-star/scale-shifts.txt}.
+It contains the three values that are necessary for modeling the star by using 
the PSF: 1. flux scaling factor, 2. x position, 3. y position.
+Take them with the following lines.
+In order to avoid the first line that is the header (starting with a #), we 
use AWK.
+
+@example
+$ scale=$(cat single-star/scale-shifts.txt | awk '!/^#/{print $1}')
+$ xcenter=$(cat single-star/scale-shifts.txt | awk '!/^#/{print $2}')
+$ ycenter=$(cat single-star/scale-shifts.txt | awk '!/^#/{print $3}')
+@end example
+
+Once we got the necessary parameters, it is possible to subtract the star, but 
this time with the shifted (and hopefully better) position.
+The last line will open the previous subtracted image and the new subtracted 
(shifted) image.
+
+@example
+$ astscript-psf-subtract label/67510-seg.fits \
+           --mode=wcs \
+           --psf=psf.fits \
+           --scale=$scale \
+           --center=$xcenter,$ycenter \
+           --output=single-star/subtracted-shifted.fits
+
+$ astscript-fits-view single-star/subtracted.fits 
single-star/subtracted-shifted.fits \
+           --ds9center=$center --ds9mode=wcs --ds9extra="-zoom 4"
+@end example
+
+Take your time and analyze the differences, do you see any?
+Which one would you say is a better subtraction?
+Difficult to say.
+We let the decission to the user, just noticing that several steps could be 
done in order to decide: re-bin the subtracted images, study if there is a 
systemmatic offset for different stars across the image, and so on.
+Note also that here we shifted the center position in steps of 1 pixel, but if 
the position of the stars you have are accurate and precisse, you could try 
using smaller shifts or even not using them at all.
+Note also that performing this center position search is time consuming, and 
when there are several stars to be modeled and subtracted, the overhead could 
be significant.
+
+Let's consider now the brightest star in the image (which is on the top edge 
of the image).
+To do it, let's take the first star from the catalog by usin the option 
@option{--head=1}.
+The necessary commands are:
+
+@example
+$ center=$(asttable outer/67510-6-10.fits --sort phot_g_mean_mag \
+                    --column=ra,dec --head=1 \
+                    | awk '@{printf "%s,%s", $1, $2@}')
+
+$ astscript-psf-scale-factor label/67510-seg.fits \
+            --mode=wcs --quiet \
+            --psf=psf.fits \
+            --center=$center \
+            --normradii=10,15 \
+            --segment=label/67510-seg.fits \
+            --output=single-star/scale.txt 
+
+$ scale=$(cat single-star/scale.txt | awk '!/^#/@{print $1@}')
+
+$ astscript-psf-subtract label/67510-seg.fits \
+           --mode=wcs \
+           --psf=psf.fits \
+           --scale=$scale \
+           --center=$center \
+           --output=single-star/subtracted.fits
+
+$ astscript-fits-view label/67510-seg.fits single-star/subtracted.fits \
+           --ds9center=$center --ds9mode=wcs --ds9extra="-zoom 4"
+@end example
+
+
+
 You will notice that there is something wrong with this ``subtraction''!
 The box of the extended PSF is clearly visible!
 The sky noise under the box is clearly larger than the rest of the noise in 
the image.