[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master e985712d: psf-scale-factor, psf-stamp and psf-
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master e985712d: psf-scale-factor, psf-stamp and psf-subtract: revered to eaab08 |
|
Date: |
Mon, 27 May 2024 14:55:11 -0400 (EDT) |
branch: master
commit e985712d4f9694133499702549593dccf7252b18
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
psf-scale-factor, psf-stamp and psf-subtract: revered to eaab08
Until now, several changes had been made in these PSF scripts from commit
fca894c03d to cad6503d44 (and later commit 115e83794b35). But these changes
caused several problems within large projects, as described in bug #65780:
https://savannah.gnu.org/bugs/?65780
With this commit, the scripts above have been revered to commit
eaab08. This was done because it will take time to find and fix the sources
of those bugs and many projects depend on these scripts. Once those
problems are fixed, the changes that were un-done with this commit can be
returned and corrected.
---
bin/script/psf-scale-factor.sh | 764 +++++++++++++++--------------------------
bin/script/psf-stamp.sh | 344 +++++++++++--------
bin/script/psf-subtract.sh | 258 ++++++--------
doc/gnuastro.texi | 278 +++------------
4 files changed, 637 insertions(+), 1007 deletions(-)
diff --git a/bin/script/psf-scale-factor.sh b/bin/script/psf-scale-factor.sh
index ca6b5ee0..029172cc 100755
--- a/bin/script/psf-scale-factor.sh
+++ b/bin/script/psf-scale-factor.sh
@@ -47,19 +47,17 @@ center=""
keeptmp=0
tmpdir=""
segment=""
-xshifts=""
-yshifts=""
normradii=""
-nobackground=0
+widthinpix=""
+nocentering=0
sigmaclip="3,0.1"
-
-
version=@VERSION@
scriptname=@SCRIPT_NAME@
+
# Output of `--usage' and `--help':
print_usage() {
cat <<EOF
@@ -99,13 +97,13 @@ $scriptname options:
-P, --psfhdu=STR HDU/extension of the PSF image.
-O, --mode=STR Coordinates mode ('wcs' or 'img').
-c, --center=FLT,FLT Center coordinates of the object.
- -n, --normradii=INT,INT Minimum and maximum radii (in pixels)
- for computing the scaling factor value.
- -S, --segment=STR Output of Segment (OBJECTS and CLUMPS).
- -s, --sigmaclip=FLT,FLT Sigma-clip multiple and tolerance.
- -x, --xshifts=FLT,FLT,FLT Min,Step,Max (in pixels) for re-centering in x.
- -y, --yshifts=FLT,FLT,FLT Min,Step,Max (in pixels) for re-centering in y.
- -b, --nobackground Do not consider a constant background component.
+ -W, --widthinpix=INT,INT Width of the stamp in pixels.
+ -n, --normradii=INT,INT Minimum and maximum radii (in pixels)
+ for computing the scaling factor value.
+ -S, --segment=STR Output of Segment (with OBJECTS and CLUMPS).
+ -s, --sigmaclip=FLT,FLT Sigma-clip multiple and tolerance.
+ -d, --nocentering Do not do the sub-pixel centering to new pix
+ grid.
Output:
-t, --tmpdir Directory to keep temporary files.
@@ -289,6 +287,9 @@ do
-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;;
+ -W|--widthinpix) widthinpix="$2";
check_v "$1" "$widthinpix"; shift;shift;;
+ -W=*|--widthinpix=*) widthinpix="${1#*=}";
check_v "$1" "$widthinpix"; shift;;
+ -W*) widthinpix=$(echo "$1" | sed -e's/-W//');
check_v "$1" "$widthinpix"; 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;;
@@ -301,15 +302,8 @@ do
-s|--sigmaclip) sigmaclip="$2";
check_v "$1" "$sigmaclip"; shift;shift;;
-s=*|--sigmaclip=*) sigmaclip="${1#*=}";
check_v "$1" "$sigmaclip"; shift;;
-s*) sigmaclip=$(echo "$1" | sed -e's/-s//');
check_v "$1" "$sigmaclip"; shift;;
-
- -x|--xshifts) xshifts="$2";
check_v "$1" "$xshifts"; shift;shift;;
- -x=*|--xshifts=*) xshifts="${1#*=}";
check_v "$1" "$xshifts"; shift;;
- -x*) xshifts=$(echo "$1" | sed -e's/-x//');
check_v "$1" "$xshifts"; shift;;
- -y|--yshifts) yshifts="$2";
check_v "$1" "$yshifts"; shift;shift;;
- -y=*|--yshifts=*) yshifts="${1#*=}";
check_v "$1" "$yshifts"; shift;;
- -y*) yshifts=$(echo "$1" | sed -e's/-y//');
check_v "$1" "$yshifts"; shift;;
- -b|--nobackground) nobackground=1; shift;;
- -b*|--nobackground=*) on_off_option_error --background -b;;
+ -d|--nocentering) nocentering=1; shift;;
+ -d*|--nocentering=*) on_off_option_error --nocentering -d;;
# Output parameters
-k|--keeptmp) keeptmp=1; shift;;
@@ -374,10 +368,9 @@ $scriptname: no center coordinates provided. You can
specify the object's centra
EOF
exit 1
else
- ncenter=$(echo $center \
- | awk 'BEGIN{FS=","} \
- {for(i=1;i<=NF;++i) c+=$i!=""} \
- END{print c}')
+ ncenter=$(echo $center | awk 'BEGIN{FS=","} \
+ {for(i=1;i<=NF;++i) c+=$i!=""} \
+ END{print c}')
if [ x$ncenter != x2 ]; then
cat <<EOF
$scriptname: '--center' (or '-c') only takes two values, but $ncenter were
given in '$center'
@@ -386,62 +379,6 @@ EOF
fi
fi
-# If x-shifts is provided.
-if [ x"$xshifts" != x ]; then
- nxshifts=$(echo $xshifts \
- | awk 'BEGIN{FS=","} \
- {for(i=1;i<=NF;++i) c+=$i!=""} \
- END{print c}')
- if [ x$nxshifts != x3 ]; then
- cat <<EOF
-$scriptname: '--xshifts' (or '-x') takes three values (xmin,xstep,xmax), but
$nxshifts were given in '$xshifts'
-EOF
- exit 1
- fi
- xstep=""
- xstep=$(echo $xshifts \
- | awk 'BEGIN{FS=","} \
- {if ($2 <= 0.0) print "xstep <= 0.0"}')
- if [ x"$xstep" != x ]; then
- cat <<EOF
-$scriptname: the 'xstep' from the option --xshifts=xmin,xstep,xmax (or '-x')
is equal or lower than zero ($xshifts), but it has to be a value above zero
-EOF
- exit 1
- fi
-else
-# If x-shifts is not provided, set it to zero with a step different of zero
-# to prevent seq chrashing
-xshifts="-0.0,0.1,+0.0"
-fi
-
-# If y-shifts is provided.
-if [ x"$yshifts" != x ]; then
- nyshifts=$(echo $yshifts \
- | awk 'BEGIN{FS=","} \
- {for(i=1;i<=NF;++i) c+=$i!=""} \
- END{print c}')
- if [ x$nyshifts != x3 ]; then
- cat <<EOF
-$scriptname: '--yshifts' (or '-x') takes three values (xmin,xstep,xmax), but
$nyshifts were given in '$yshifts'
-EOF
- exit 1
- fi
- xstep=""
- xstep=$(echo $yshifts \
- | awk 'BEGIN{FS=","} \
- {if ($2 <= 0.0) print "ystep <= 0.0"}')
- if [ x"$xstep" != x ]; then
- cat <<EOF
-$scriptname: the 'ystep' from the option --yshifts=xmin,ystep,xmax (or '-y')
is equal or lower than zero ($yshifts), but it has to be a value above zero
-EOF
- exit 1
- fi
-else
-# If y-shifts is not provided, set it to zero with a step different of zero
-# to prevent seq chrashing
-yshifts="-0.0,0.1,+0.0"
-fi
-
# If a normalization range is not given at all.
if [ x"$normradii" = x ]; then
cat <<EOF
@@ -529,20 +466,6 @@ fi
-# Compute the width of the stamp
-# ------------------------------
-#
-# If the width of the stamp is specified by the user, use that size.
-# Otherwise, make the stamp width of the same size than two times the
-# maximum radius for computing the flux factor. This is the maximum radius
-# that is needed for computing the scaling flux factor.
-xwidth=$(astarithmetic $normradiusmax float32 2.0 x 1.0 + --quiet)
-ywidth=$(astarithmetic $normradiusmax float32 2.0 x 1.0 + --quiet)
-
-
-
-
-
# Transform WCS to IMG center coordinates
# ---------------------------------------
#
@@ -551,449 +474,314 @@ ywidth=$(astarithmetic $normradiusmax float32 2.0 x 1.0
+ --quiet)
# 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)
+ xcenter=$(echo "$xycenter" | awk '{print $1}')
+ ycenter=$(echo "$xycenter" | awk '{print $2}')
+else
+ xcenter=$xcoord
+ ycenter=$ycoord
+fi
- 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. The fake WCS is constructed as
- # follows: 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
- fake_wcs=$tmpdir/fake-wcs-$objectid.fits
- 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}')
- 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}')
-psfcenter2=$(echo $psfnaxis2 | awk '{print $1/2+0.5}')
-
-# Create the fake WCS image. No data, just WCS.
-psf_fake_wcs=$tmpdir/psf-fake-wcs.fits
-echo "1 $psfcenter1 $psfcenter2 1 1 1 1 1 1 1" \
- | astmkprof --oversample=1 \
- --type=uint8 \
- --output=$psf_fake_wcs \
- --cdelt=$cdelt1,$cdelt2 \
- --crval=$crval1,$crval2 \
- --crpix=$psfcenter1,$psfcenter2 \
- --mergedsize=$psfnaxis1,$psfnaxis2
-
-# Inject the WCS header into the image.
-psf_wcs=$tmpdir/psf-wcsed.fits
-astarithmetic $psf --hdu=$psfhdu \
- --wcsfile=$psf_fake_wcs --output=$psf_wcs
-
-# Transform the center in img coordinates to center in wcs
-psfxywcs=$(echo "$psfcenter1,$psfcenter2" \
- | asttable --column='arith $1 $2 img-to-wcs' \
- --wcsfile=$psf_wcs --wcshdu=1 $quiet \
- | awk '{printf "%.10e,%.10e\n", $1, $2}')
-
-# Warp the PSF
-# ------------
+# Compute the width of the stamp
+# ------------------------------
#
-# 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
+# If the width of the stamp is specified by the user, use that size.
+# Otherwise, make the stamp width of the same size than two times the
+# maximum radius for computing the flux factor. This is the maximum radius
+# that is needed for computing the flux value.
+if [ x"$widthinpix" = x ]; then
+ xwidthinpix=$(astarithmetic $normradiusmax float32 2.0 x 1.0 + --quiet)
+ ywidthinpix=$(astarithmetic $normradiusmax float32 2.0 x 1.0 + --quiet)
+else
+ xwidthinpix=$(echo $widthinpix | awk '{print $1}')
+ ywidthinpix=$(echo $widthinpix | awk '{print $2}')
+fi
-# Build a radial profile image
-# ----------------------------
+# Crop the original image around the object
+# -----------------------------------------
#
-# It will be used to only select pixels within the requested radial range.
-xradcenter=$(astfits $psf_warped --keyval NAXIS1 --quiet \
- | awk '{print $1/2+0.5}')
-yradcenter=$(astfits $psf_warped --keyval NAXIS2 --quiet \
- | awk '{print $1/2+0.5}')
-maxradius=$(printf "$xwidth\n$ywidth" \
- | aststatistics --maximum --quiet)
-rad_warped=$tmpdir/warped-radial-$objectid.fits
-echo "1 $xradcenter $yradcenter 7 $maxradius 0 0 1 1 1" \
- | astmkprof --background=$psf_warped --clearcanvas \
- --oversample=1 --output=$rad_warped $quiet
+# Crop the object around its center with the given stamp size width.
+cropped=$tmpdir/cropped-$objectid.fits
+astcrop $inputs --hdu=$hdu --mode=img \
+ --center=$xcenter,$ycenter \
+ --width=$xwidthinpix,$ywidthinpix\
+ --output=$cropped $quiet
+
+# If the cropped image is not generated, it may happen that it does not
+# overlap with the input image. Save a nan value as the output and
+# continue.
+if ! [ -f $cropped ]; then
+
+ multifactor=nan
+
+ # Let the user know what happened.
+ if [ x"$quiet" = x ]; then all_nan_warning; fi
+
+ # Print the multiplication factor in standard output or in a given file
+ if [ x"$output" = x ]; then echo $multifactor
+ else echo $multifactor > $output
+ fi
+ # If the user does not specify to keep the temporal files with the option
+ # `--keeptmp', then remove the whole directory.
+ if [ $keeptmp = 0 ]; then
+ rm -r $tmpdir
+ fi
+ exit 0
+fi
-# Shifts grid
-# -----------
+# Crop and unlabel the segmentation image
+# ---------------------------------------
#
-# Here, a grid of center positions are defined. The image will be warped
-# using each center. At the end, the one that has the smaller standard
-# deviation within the ring of computing the flux scaling factor will be
-# provided as the best option. The shifts are specified in pixels, so, they
-# need to be transformed to WCS shifts by using the pixel scale in deg/pix.
-pscales=$(astfits $inputs_wcs --hdu $hdu --pixelscale --quiet)
-xpscale=$(echo $pscales | awk '{print $1}')
-ypscale=$(echo $pscales | awk '{print $2}')
-xshiftswcs=$(echo $xshifts \
- | awk -v p="$xpscale" -F ","\
- '{printf "%.10f %.10f %.10f\n", p*$1,p*$2,p*$3}')
-yshiftswcs=$(echo $yshifts \
- | awk -v p="$ypscale" -F ","\
- '{printf "%.10f %.10f %.10f\n", p*$1,p*$2,p*$3}')
-
-# For each offset in the x and y direction
-for xs in $(seq $xshiftswcs); do
- for ys in $(seq $yshiftswcs); do
-
- # Compute the shifts in WCS
- xspix=$(astarithmetic $xs $xpscale / --quiet \
- | awk '{printf "%.10f\n", $1}')
- yspix=$(astarithmetic $ys $ypscale / --quiet \
- | awk '{printf "%.10f\n", $1}')
-
- # New shifted center in WCS
- xcwcs_shift=$(astarithmetic $xcwcs $xs + --quiet)
- ycwcs_shift=$(astarithmetic $ycwcs $ys + --quiet)
-
- #To ease the filename reading, put the offsets in pixels
- label_shift=xs_ys_"$xspix"_"$yspix"
-
-
- # Warp (WCS) the original image around the object
- # -----------------------------------------------
- #
- # Warp the object around its center with the given stamp size width.
- warped=$tmpdir/warped-$objectid-$label_shift.fits
- astwarp $inputs_wcs \
- --hdu=$hdu\
- --widthinpix \
- --width=$xwidth,$ywidth \
- --output=$warped $quiet \
- --center=$xcwcs_shift,$ycwcs_shift
-
-
- # If the warped image is not generated, it may happen that it does not
- # overlap with the input image. Save a nan value as the output and
- # continue.
- if ! [ -f $warped ]; then
-
- outvalues="nan nan nan"
-
- # Let the user know what happened.
- if [ x"$quiet" = x ]; then all_nan_warning; fi
-
- # Print the values in standard output or in a given file
- if [ x"$output" = x ]; then echo $outvalues
- else
- echo "# factor xcenter ycenter" > $output
- echo $outvalues >> $output
-
- fi
-
- # If the user does not specify to keep the temporal files with the
option
- # `--keeptmp', then remove the whole directory.
- if [ $keeptmp = 0 ]; then
- rm -r $tmpdir
- fi
-
- exit 0
-
- fi
-
-
- # Crop and unlabel the segmentation image
- # ---------------------------------------
- #
- # If the user provides a segmentation image, treat it appropiately
- # in order to mask all objects that are not the central one. If
- # not, just consider that the cropped and masked image is the
- # cropped (not masked) image. The process is as follow:
- # - Compute the central clump and object labels. This is done
- # with the option '--oneelemstdout' of Crop.
- # - 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.
- # - In the original cropped mask, convert pixels belonging to the
- # 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"$segment" != x ]; then
-
- # Find the object and clump labels of the target.
- clab=$(astcrop $segment --hdu=CLUMPS --mode=img --width=1,1 \
- --oneelemstdout --center=$xcimg,$ycimg --quiet)
- olab=$(astcrop $segment --hdu=OBJECTS --mode=img --width=1,1 \
- --oneelemstdout --center=$xcimg,$ycimg --quiet)
-
- # If for any reason, a clump or object label couldn't be
- # initialized at the given coordiante, simply ignore this
- # step. But print a warning so the user is informed of the
- # situation (and that this is a bug: 'clab' should be
- # initialized!).
- if [ x"$clab" = x -o x"$olab" = x ]; then
- if [ x"$quiet" = x ]; then
- cat <<EOF
+# If the user provides a segmentation image, treat it appropiately in order
+# to mask all objects that are not the central one. If not, just consider
+# that the cropped and masked image is the cropped (not masked) image. The
+# process is as follow:
+# - Compute the central clump and object labels. This is done with the
+# option '--oneelemstdout' of Crop.
+# - 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 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.
+# - Mask all non zero pixels in the mask image as nan values.
+if [ x"$segment" != x ]; then
+
+ # Find the object and clump labels of the target.
+ clab=$(astcrop $segment --hdu=CLUMPS --mode=img --width=1,1 \
+ --oneelemstdout --center=$xcenter,$ycenter \
+ --quiet)
+ olab=$(astcrop $segment --hdu=OBJECTS --mode=img --width=1,1 \
+ --oneelemstdout --center=$xcenter,$ycenter \
+ --quiet)
+
+ # If for any reason, a clump or object label couldn't be initialized at
+ # the given coordiante, simply ignore this step. But print a warning so
+ # the user is informed of the situation (and that this is a bug: 'clab'
+ # should be initialized!).
+ if [ x"$clab" = x -o x"$olab" = x ]; then
+ if [ x"$quiet" = x ]; then
+ cat <<EOF
$scriptname: WARNING: no clump or object label could be found in '$segment'
for the coordinate $center
EOF
- fi
- warped_masked=$warped
- else
- # To help in debugging (when '--quiet' is not called)
- if [ x"$quiet" = x ]; then
- echo "$scriptname: $segment: at $center, found clump $clab
in object $olab"
- fi
-
- # In order to warp the objects and clumps image, WCS
- # information is needed. If mode is wcs, just copy the HDU.
- # If mode is img (no wcs info), it is necessary to inject
- # the fake WCS information in advance.
- clumps_wcs=$tmpdir/clumps-wcs-$objectid-$label_shift.fits
- objects_wcs=$tmpdir/objects-wcs-$objectid-$label_shift.fits
- if [ "$mode" = wcs ]; then
- astfits $segment --copy CLUMPS --output $clumps_wcs
- astfits $segment --copy OBJECTS --output $objects_wcs
- else
- astarithmetic $segment --hdu CLUMPS \
- --wcsfile=$fake_wcs --output=$clumps_wcs
- astarithmetic $segment --hdu OBJECTS \
- --wcsfile=$fake_wcs --output=$objects_wcs
- fi
-
- # 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=$croped_clp $quiet \
- --center=$xcwcs_shift,$ycwcs_shift
- astcrop $objects_wcs \
- --widthinpix \
- --width=$xwidth,$ywidth \
- --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 \
- $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 \
- i omask cmask or nan where \
- --output=$warped_masked $quiet
- fi
- else
- warped_masked=$warped
fi
-
- # Set the used the center position for the output. If the original
- # mode requested was IMG, the current center postion (that is in
- # WCS) needs to be transformed to IMG.
- xc_out=$xcwcs_shift
- yc_out=$ycwcs_shift
- if [ x"$mode" = x"img" ]; then
- xcyc_shift=$(echo "$xcwcs_shift,$ycwcs_shift" \
- | asttable --column='arith $1 $2 wcs-to-img' \
- --wcsfile=$inputs_wcs --wcshdu=1 $quiet)
- xc_out=$(echo "$xcyc_shift" | awk '{print $1}')
- yc_out=$(echo "$xcyc_shift" | awk '{print $2}')
+ cropped_masked=$cropped
+ else
+ # To help in debugging (when '--quiet' is not called)
+ if [ x"$quiet" = x ]; then
+ echo "$scriptname: $segment: at $center, found clump $clab in
object $olab"
fi
- # Mask all but the wanted pixels of the ring.
- multipimg=$tmpdir/for-factor-$objectid-$label_shift.fits
- astarithmetic $warped_masked -h1 set-i \
- $psf_warped -h1 set-p \
- $rad_warped -h1 set-r \
- r $normradiusmin lt r $normradiusmax ge or set-m \
- i p / m nan where --output $multipimg $quiet
-
- # Find the multiplication factor, and the STD
- stats=$(aststatistics $multipimg --quiet \
- --sclipparams=$sigmaclip \
- --sigclip-median --sigclip-std)
-
- # Collect the data: center position, flux scale, and shifts.
- values=$tmpdir/stats-$objectid-$label_shift.txt
- echo "#xcenter ycenter fscale fscalestd xshift yshift xshiftpix
yshiftpix" > $values
- echo "$xc_out $yc_out $stats $xs $ys $xspix $yspix
" >> $values
-
- done
-done
+ # 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 \
+ --width=$xwidthinpix,$ywidthinpix \
+ --center=$xcenter,$ycenter \
+ --output=$cropobj $quiet
+ astcrop $segment --hdu=CLUMPS --mode=img \
+ --width=$xwidthinpix,$ywidthinpix \
+ --center=$xcenter,$ycenter \
+ --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 $quiet
+ fi
+else
+ cropped_masked=$cropped
+fi
-# Find the shifted position that gives the smallest STD value
-# -----------------------------------------------------------
+# Function to simplify checks in each dimension
+# ---------------------------------------------
#
-# Once all data has been computed for all shifted positions, here the one
-# that gives the smallest STD value is obtained. This is the one that will
-# be given as the best output.
+# If the bottom-left edge of the crop is outside the image, the minimum
+# coordinate needs to be modified from the value in the 'ICF1PIX' keyword
+# of the output of the Crop program.
+first_pix_in_img () {
+
+ # This function takes a single argument: the dimension to work on.
+ dim=$1
+
+ # Find the minimim and maximum from the cropped image.
+ if [ $dim = x ]; then
+ width=$xwidthinpix
+ min=$(echo "$overlaprange" | awk '{print $1}')
+ max=$(echo "$overlaprange" | awk '{print $2}')
+ elif [ $dim = y ]; then
+ width=$ywidthinpix
+ min=$(echo "$overlaprange" | awk '{print $3}')
+ max=$(echo "$overlaprange" | awk '{print $4}')
+ else
+ cat <<EOF
+$scriptname: a bug! Please contact with us at bug-gnuastro@gnu.org. The value
of 'dim' is "$dim" and not recognized
+EOF
+ exit 1
+ fi
-# Stack data of all shifted images and sort them by the STD (4th) column
-stats_all=$tmpdir/stats-all-$objectid.fits
-cat $tmpdir/stats-$objectid-*.txt \
- | asttable --sort=4 --output $stats_all
+ # Check if it is necessary to correct the position. If the bottom-left
+ # corner of the image (which defines the zero coordinate along both
+ # dimensions) has fallen outside of the image, then the minimum value
+ # of 'ICF1PIX' in that dimension will be '1'.
+ #
+ # Therefore, if the minimum value is not '1', then everything is fine
+ # and no correction is necessary! We'll just return the minimum
+ # value. Otherwise, we need to return a negative minimum pixel value to
+ # obtain the correct position of the object from 'ICF1PIX' (which only
+ # has the region of the input that overlapped with the output). It is
+ # necessary to add '1' because the first pixel has a coordiante of 1,
+ # not 0.
+ if [ $min = 1 ]; then
+ echo "$max $width" | awk '{if($1==$2) print 1; else print $1-$2+1}'
+ else
+ echo $min
+ fi
+}
-# Get the shift with the minimum STD value
-minstd=$(asttable $stats_all -c4 --head=1 --quiet)
-# Get the coordinates
-xcyc_best=$(asttable $stats_all --equal=4,$minstd -c1,2 --quiet)
+# In the FITS standard, the integer coordinates are defined on the center
+# of each pixel. On the other hand, the centers of objects can be anywhere
+# (not exactly in the center of the pixel!). In this scenario, we should
+# move the center of the object to the center of the pixel with a sub-pixel
+# warp. In this part of the script we will calculate the necessary
+# sub-pixel translation and use Warp to center the object on the pixel
+# grid. The user can disable this with the '--nocentering' option.
+if [ $nocentering = 0 ]; then
+
+ # Read the overlap range from the 'ICF1PIX' keyword (which is printed
+ # in all outputs of Crop).
+ overlaprange=$(astfits $cropped -h0 --keyvalue=ICF1PIX --quiet \
+ | sed -e's|:| |g' -e's|,| |')
+
+ # Calculate the position of the bottom-left pixel of the cropped image
+ # in relation to the input image.
+ minx=$(first_pix_in_img x)
+ miny=$(first_pix_in_img y)
+
+ # Due to cropping, the image coordinate should be shifted by a certain
+ # number (integer) of pixels (leaving the sub-pixel center intact).
+ stmcoord=$(echo "$xcenter $ycenter $minx $miny" \
+ | awk '{printf "%g %g\n", $1-$3+1, $2-$4+1}')
+
+ # Calculate the displacement (or value to give to '--translate' in
+ # Warp), and the new center of the star after translation (or value of
+ # '--center' in Crop):
+ # 1. At first we extract the sub-pixel displacement ('x' and 'y').
+ # 2. If the displacement is larger than 0.5, then after the warp, the
+ # center will be shifted by one pixel. Otherwise, if the
+ # displacement is smaller than 0.5, we shift the grid backwards by
+ # mutiplying it by -1 and don't add any shift.
+ # 3. Add 'xshift' and 'yshit' valut to integer of x and y coordinate
+ # and add final value to one.
+ warpcoord=$(echo "$stmcoord" \
+ | awk '{x=$1-int($1); y=$2-int($2); \
+ if(x>0.5){x=1-x; xshift=1} else {x*=-1; xshift=0}; \
+ if(y>0.5){y=1-y; yshift=1} else {y*=-1; yshift=0}; \
+ printf("%f,%f %d,%d\n", x, y, \
+ int($1)+xshift+1, int($2)+yshift+1)}')
+
+ # Warp image based on the measured displacement (first component of the
+ # output above).
+ warpped=$tmpdir/cropped-masked-warpforcenter-$objectid.fits
+ DXY=$(echo "$warpcoord" | awk '{print $1}')
+ astwarp $cropped_masked --translate=$DXY --output=$warpped $quiet
+
+ # Crop image based on the calculated shift (second component of the
+ # output above).
+ centermsk=$tmpdir/cropped-masked-centered-$objectid.fits
+ CXY=$(echo "$warpcoord" | awk '{print $2}')
+ astcrop $warpped -h1 \
+ --mode=img --output=$centermsk $quiet \
+ --center=$CXY --width=$xwidthinpix,$ywidthinpix
+else
+ # If the user did not want to correct the center of image, we'll use
+ # the raw crop as input for the next step.
+ centermsk=$cropped_masked
+fi
-# Get the flux scaling factor (after background estimation)
-# ---------------------------------------------------------
-#
-# Up to here, the coordinates that give the smallest STD have been
-# computed. Still it is necessary to compute the flux factor. There hare
-# two options now for computing it.
-# - If a background component is NOT considered, then just read the flux
-# factor that is already computed.
-# - If a background component IS considered, estimate it by using the
-# radial profiles, subtract this value, and estimate the flux factor.
-#
-if [ x$nobackground = x1 ]; then
- back_val=0.0
- factor_val=$(asttable $stats_all --equal=4,$minstd -c3 --quiet)
-else
- # Construct the label to recover the best warped image
- label_shift_best=$(asttable $stats_all --equal=4,$minstd -c7,8 --quiet \
- | awk '{printf "%.10f_%.10f\n", $1,$2}')
-
- # This is the warped and shifted image with the best (lowest) STD
- warped_best=$tmpdir/warped-$objectid-xs_ys_$label_shift_best.fits
-
- # Generate the radial profiles
- rprofile_psf=$tmpdir/rprofile-psf-$objectid-$label_shift_best.txt
- rprofile_img=$tmpdir/rprofile-img-$objectid-$label_shift_best.txt
- astscript-radial-profile $psf_warped --output $rprofile_psf
- astscript-radial-profile $warped_best --output $rprofile_img
-
- # Read the number of lines, and the contents of the tables into arrays,
- # skipping lines starting with #
- n=$(asttable $rprofile_img | wc -l)
- IFS=$'\n' array_img=($(grep -v '^#' $rprofile_img))
- IFS=$'\n' array_psf=($(grep -v '^#' $rprofile_psf))
-
- mixed_table=$tmpdir/mixed-star-psf-$objectid-$label_shift_best.fits
- for i in $(seq 0 $(($n - 1))); do
- for j in $(seq 0 $(($n - 1))); do
- printf "%s %s %s %s %s %s\n" ${array_img[i]} ${array_img[j]}
${array_psf[i]} ${array_psf[j]}
- done
- done | asttable -c1,3,2,4,6,8 --output $mixed_table
-
- for_back_fits=$tmpdir/for-background-stats.fits
- asttable $mixed_table \
- -c1,2,3,4,5,6 \
- -c'arith $3 $4 - $5 $6 - /' \
- -c'arith $3 $5 $3 $4 - $5 $6 - / x -' \
- --output $for_back_fits
-
- # Find the background value
- back_val=$(aststatistics $for_back_fits -c8 --quiet \
- --sclipparams=$sigmaclip --sigclip-mean)
-
- # Mask all but the wanted pixels of the ring. Subtract the
- # background value computed above (or vback=0.0 if not estimated)
- multipimg=$tmpdir/for-factor-$objectid-$label_shift_best.fits
- astarithmetic $warped_masked -h1 $back_val - set-i \
- $psf_warped -h1 set-p \
- $rad_warped -h1 set-r \
- r $normradiusmin lt r $normradiusmax ge or set-m \
- i p / m nan where --output $multipimg $quiet
-
- # Find the multiplication factor and its STD
- factor_val=$(aststatistics $multipimg --quiet \
- --sclipparams=$sigmaclip --sigclip-median)
-fi
+# Crop the PSF image with the same size.
+psfcropped=$tmpdir/cropped-psf-$objectid.fits
+psfxcenter=$(astfits $psf -h$psfhdu --keyvalue=NAXIS1 --quiet \
+ | awk '{print $1/2+1}')
+psfycenter=$(astfits $psf -h$psfhdu --keyvalue=NAXIS2 --quiet \
+ | awk '{print $1/2+1}')
+astcrop $psf --hdu=$psfhdu --mode=img \
+ --center=$psfxcenter,$psfycenter \
+ --width=$xwidthinpix,$ywidthinpix \
+ --output=$psfcropped $quiet
+
+
+
+
+
+# Build a radial profile image. It will be used to only select pixels
+# within the requested radial range.
+xradcenter=$(echo $xwidthinpix | awk '{print $1/2+1}')
+yradcenter=$(echo $ywidthinpix | awk '{print $1/2+1}')
+maxradius=$(printf "$xwidthinpix\n$ywidthinpix" \
+ | aststatistics --maximum --quiet)
+radcropped=$tmpdir/cropped-radial-$objectid.fits
+echo "1 $xradcenter $yradcenter 7 $maxradius 0 0 1 1 1" \
+ | astmkprof --background=$psfcropped --clearcanvas \
+ --oversample=1 --output=$radcropped $quiet
+
+
+
+
+
+# Find the multiplication factor
+multipimg=$tmpdir/for-factor-$objectid.fits
+astarithmetic $centermsk -h1 set-i \
+ $psfcropped -h1 set-p \
+ $radcropped -h1 set-r \
+ r $normradiusmin lt r $normradiusmax ge or set-m \
+ i p / m nan where --output $multipimg $quiet
+multifactor=$(aststatistics $multipimg --sigclip-median \
+ --sclipparams=$sigmaclip --quiet)
+
-# Output
-# ------
-#
-# Collect the output data and print the (standard output or write a file)
-outvalues="$factor_val $xcyc_best"
-if [ x"$output" = x ]; then
- echo $outvalues
-else
- echo "# factor xcenter ycenter" > $output
- echo $outvalues >> $output
- echo; echo "Output written to '$output'."
+# Print the multiplication factor in standard output or in a given file
+if [ x"$output" = x ]; then echo $multifactor
+else echo $multifactor > $output
fi
diff --git a/bin/script/psf-stamp.sh b/bin/script/psf-stamp.sh
index 25dc5d86..b1138326 100644
--- a/bin/script/psf-stamp.sh
+++ b/bin/script/psf-stamp.sh
@@ -1,7 +1,7 @@
#!/bin/sh
# Construct a star stamp to build a Point Spread Function (PSF). This
-# script will consider a center position and then it will warp the original
+# script will consider a center position and then it will crop the original
# image around that center with the specified size. After that, it will
# compute the radial profile to obtain a normalization value. Finally, the
# output stamp will be normalized by dividing by the normalization value.
@@ -60,6 +60,7 @@ axisratio=1
normradii=""
sigmaclip=""
widthinpix=""
+nocentering=0
normop="median"
positionangle=0
version=@VERSION@
@@ -114,6 +115,7 @@ $scriptname options:
-Q, --axis-ratio=FLT Axis ratio for ellipse maskprofile (A/B).
-p, --position-angle=FLT Position angle for ellipse mask profile.
-s, --sigmaclip=FLT,FLT Sigma-clip multiple and tolerance.
+ -d, --nocentering Don't warp the center coord to central pixel.
Output:
-o, --output Output table with the radial profile.
@@ -315,9 +317,9 @@ do
-Q|--axis-ratio) axisratio="$2";
check_v "$1" "$axisratio"; shift;shift;;
-Q=*|--axis-ratio=*) axisratio="${1#*=}";
check_v "$1" "$axisratio"; shift;;
-Q*) axisratio=$(echo "$1" | sed -e's/-Q//');
check_v "$1" "$axisratio"; shift;;
- -p|--position-angle) positionangle="$2";
check_v "$1" "$positionangle"; shift;shift;;
- -p=*|--position-angle=*) positionangle="${1#*=}";
check_v "$1" "$positionangle"; shift;;
- -p*) positionangle=$(echo "$1" | sed -e's/-p//');
check_v "$1" "$positionangle"; shift;;
+ -p|--position-angle) positionangle="$2";
check_v "$1" "$positionangle"; shift;shift;;
+ -p=*|--position-angle=*) positionangle="${1#*=}";
check_v "$1" "$positionangle"; shift;;
+ -p*) positionangle=$(echo "$1" | sed
-e's/-p//'); check_v "$1" "$positionangle"; shift;;
# Output parameters
-k|--keeptmp) keeptmp=1; shift;;
@@ -325,6 +327,8 @@ do
-t|--tmpdir) tmpdir="$2"; check_v "$1"
"$tmpdir"; shift;shift;;
-t=*|--tmpdir=*) tmpdir="${1#*=}"; check_v "$1"
"$tmpdir"; shift;;
-t*) tmpdir=$(echo "$1" | sed -e's/-t//'); check_v "$1"
"$tmpdir"; shift;;
+ -d|--nocentering) nocentering=1; shift;;
+ -d*|--nocentering=*) on_off_option_error --nocentering -d;;
-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;;
@@ -533,86 +537,46 @@ fi
-# Compute the center in WCS for Warping
-# -------------------------------------
+# Transform WCS to IMG center coordinates
+# ---------------------------------------
#
# If the original coordinates have been given in WCS or celestial units
-# (RA/DEC) they will be used by Warp. Otherwise, here fake WCS information
-# is generated to be able to use Warp in the next steps. This is necessary
-# for allocating the object with sub-pixel precission at the center of the
-# output image. The IMG coordinates of the center are also computed because
-# in next steps they are necessary if a segment image is provided (to
-# compute the object and clumps label).
+# (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
-
+ xycenter=$(echo "$xcoord,$ycoord" \
+ | asttable --column='arith $1 $2 wcs-to-img' \
+ --wcsfile=$inputs --wcshdu=$hdu $quiet)
+ xcenter=$(echo "$xycenter" | awk '{print $1}')
+ ycenter=$(echo "$xycenter" | awk '{print $2}')
else
- # Fake WCS information to be injected as header to the original image.
- # The center is assumed to be at RA=180,DEC=0 with a pixel scale of 1
- # arcsec in both axeses. Warp will interpret this for warping.
- crval1=180.0
- crval2=0.0
- cdelt1=1.0
- cdelt2=1.0
- naxis1=$(astfits $inputs --hdu=$hdu --keyval NAXIS1 --quiet)
- naxis2=$(astfits $inputs --hdu=$hdu --keyval NAXIS2 --quiet)
- fake_wcs=$tmpdir/fake-wcs-$objectid.fits
- echo "1 $xcoord $ycoord 1 1 1 1 1 1 1" \
- | astmkprof --type=uint8 \
- --oversample=1 \
- --output=$fake_wcs \
- --cunit="arcsec,arcsec" \
- --cdelt=$cdelt1,$cdelt2 \
- --crval=$crval1,$crval2 \
- --crpix=$xcoord,$ycoord \
- --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}')
-
+ xcenter=$xcoord
+ ycenter=$ycoord
fi
-# Warp the original image around the object
+# Crop the original image around the object
# -----------------------------------------
#
-# Warp the object around its center with the given stamp size width. It may
+# Crop the object around its center with the given stamp size width. It may
# happen that the given coordinate is fully outside of the image (within
-# the requested stamp-width). In this case Warp won't generate any output,
-# so we are checking the existance of the '$warped' file. If not created
+# the requested stamp-width). In this case Crop won't generate any output,
+# so we are checking the existance of the '$cropped' file. If not created
# we will create a fully NaN-valued image that can be ignored in stacks.
-warped=$tmpdir/warped-$objectid.fits
-astwarp $inputs_wcs \
- --hdu=1 \
- --widthinpix \
- --center=$xcwcs,$ycwcs \
- --output=$warped $quiet \
- --width=$xwidthinpix,$ywidthinpix
-if ! [ -f $warped ]; then
+#
+# We are setting '--checkcenter=0' in Crop so it doesn't check if the
+# central pixel is covered or not (we may be interested in the outer parts
+# of the PSF of a star that is centered outside of the image).
+cropped=$tmpdir/cropped.fits
+astcrop $inputs --hdu=$hdu --mode=img --checkcenter=0 \
+ --width=$xwidthinpix,$ywidthinpix \
+ --center=$xcenter,$ycenter \
+ --output=$cropped $quiet
+if ! [ -f $cropped ]; then
# 'makenew' will generate a fully 0-valued image. Adding (with '+') any
# number by NaN will be NaN. Therefore the output of the Arithmetic
@@ -634,83 +598,70 @@ fi
-
-# Warp and unlabel the segmentation image
+# Crop and unlabel the segmentation image
# ---------------------------------------
#
# If the user provides a segmentation image, treat it appropiately in order
# to mask all objects that are not the central one. If not, just consider
-# that the warped and masked image is the warped (not masked) image. The
+# that the cropped and masked image is the cropped (not masked) image. The
# process is as follow:
# - Compute the central clump and object labels. This is done with the
# option '--oneelemstdout' of Crop.
-# - Warp the CLUMPS image.
-# - Warp the OBJECTS image.
-# - On the original warped image, convert pixels belonging to the central
-# object to zeros. By doing this, the central object becomes as sky.
+# - 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 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.
# - Mask all non zero pixels in the mask image as nan values.
if [ x"$segment" != x ]; then
+
# Find the object and clump labels of the target.
clab=$(astcrop $segment --hdu=CLUMPS --mode=img --width=1,1 \
- --oneelemstdout --center=$xcimg,$ycimg --quiet)
+ --oneelemstdout --center=$xcenter,$ycenter \
+ --quiet)
olab=$(astcrop $segment --hdu=OBJECTS --mode=img --width=1,1 \
- --oneelemstdout --center=$xcimg,$ycimg --quiet)
+ --oneelemstdout --center=$xcenter,$ycenter \
+ --quiet)
# If for any reason, a clump or object label couldn't be initialized at
# the given coordiante, simply ignore this step. But print a warning so
# the user is informed of the situation (and that this is a bug: 'clab'
# should be initialized!).
if [ x"$clab" = x -o x"$olab" = x ]; then
- if [ x"$quiet" = x ]; then
- cat <<EOF
+ cat <<EOF
$scriptname: WARNING: no clump or object label could be found in '$segment'
for the coordinate $center
EOF
- fi
- warped_masked=$warped
+ cropped_masked=$cropped
+
else
# To help in debugging (when '--quiet' is not called)
- if [ x"$quiet" = x ]; then
+ if [ x"$quiet" = x ] && [ x"$center" != x ]; then
echo "$scriptname: $segment: at $center, found clump $clab in
object $olab"
fi
- # In order to warp the objects and clumps image, WCS
- # information is needed. If mode is wcs, just copy the HDU.
- # If mode is img (no wcs info), it is necessary to inject
- # the fake WCS information in advance.
- clumps_wcs=$tmpdir/clumps-wcs-$objectid.fits
- objects_wcs=$tmpdir/objects-wcs-$objectid.fits
- if [ "$mode" = wcs ]; then
- astfits $segment --copy CLUMPS --output $clumps_wcs
- astfits $segment --copy OBJECTS --output $objects_wcs
- else
- astarithmetic $segment --hdu CLUMPS \
- --wcsfile=$fake_wcs --output=$clumps_wcs
- astarithmetic $segment --hdu OBJECTS \
- --wcsfile=$fake_wcs --output=$objects_wcs
- fi
-
- # Maske the detection map of Objects and Clumps of target in image.
+ # Crop the object and clump image to same size as desired stamp.
+ cropclp=$tmpdir/cropped-clumps.fits
+ cropobj=$tmpdir/cropped-objects.fits
+ astcrop $segment --hdu=OBJECTS --mode=img \
+ --center=$xcenter,$ycenter \
+ --width=$xwidthinpix,$ywidthinpix --output=$cropobj $quiet
+ astcrop $segment --hdu=CLUMPS --mode=img \
+ --center=$xcenter,$ycenter \
+ --width=$xwidthinpix,$ywidthinpix --output=$cropclp $quiet
+
+ # Mask all the undesired regions.
cropped_masked=$tmpdir/cropped-masked.fits
- astarithmetic $segment --hdu=INPUT-NO-SKY set-i \
- $segment --hdu=CLUMPS set-c \
- $segment --hdu=OBJECTS set-o \
+ astarithmetic $cropped --hdu=1 set-i --output=$cropped_masked $quiet \
+ $cropobj --hdu=1 set-o \
+ $cropclp --hdu=1 set-c \
\
c o $olab ne 0 where c $clab eq -1 where 0 gt \
- 1 dilate set-mask \
+ 1 dilate set-cmask \
o o $olab eq 0 where set-omask \
- mask omask mask or nan where -o$cropped_masked
-
- # Warp the mask image.
- cropped_masked_warp=$tmpdir/cropped-masked-warp.fits
- astwarp $cropped_masked \
- --gridfile=$warped $quiet \
- --output=$cropped_masked_warp
-
- # Mask extra objects in warp image.
- warped_masked=$tmpdir/warped-masked-$objectid.fits
- astarithmetic $warped -h1 set-w \
- $cropped_masked_warp -h1 set-m \
- w m isblank nan where -o$warped_masked
+ i omask cmask or nan where
fi
# Apply the signal-to-noise threshold if the user provided one, and
@@ -720,48 +671,163 @@ EOF
# Mask all the pixels that are below the threshold based on the
# signal to noise value which is obtained form the SKY_STD. Finally
# label the separate regions.
- warpsnt=$tmpdir/warped-snt.fits
- warpstd=$tmpdir/warped-std.fits
- warplab=$tmpdir/warped-lab.fits
- astwarp $segment --hdu=SKY_STD \
- --gridfile=$warped \
- --output=$warpstd $quiet
+ cropsnt=$tmpdir/cropped-snt.fits
+ cropstd=$tmpdir/cropped-std.fits
+ croplab=$tmpdir/cropped-lab.fits
+ astcrop $segment --hdu=SKY_STD --mode=img \
+ --center=$xcenter,$ycenter \
+ --width=$xwidthinpix,$ywidthinpix \
+ --output=$cropstd $quiet
# Fill the NAN pixels with maximum value of the image plus one (so
# the fill value is not present in the image).
- fillval=$(astarithmetic $warped_masked maxvalue 1 + -q)
- fill=$tmpdir/warped-masked-fill-nan-pix.fits
- astarithmetic $warped_masked set-i i i isblank $fillval \
+ fillval=$(astarithmetic $cropped_masked maxvalue 1 + -q)
+ fill=$tmpdir/cropped-masked-fill-nan-pix.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 \
- $warpstd -h1 set-s \
- v s / set-sn \
+ astarithmetic $fill -h1 set-v \
+ $cropstd -h1 set-s \
+ v s / set-sn \
v sn $snthresh lt \
2 dilate nan where set-all \
- all tofilefree-$warpsnt \
+ all tofilefree-$cropsnt \
all isnotblank 2 connected-components \
- --output=$warplab
+ --output=$croplab
# Extract the label of the central coordinate.
- id=$(astcrop $warplab -h1 --mode=$mode \
- --center=$xcimg,$ycimg \
+ id=$(astcrop $croplab -h1 --mode=$mode \
+ --center=$xcoord,$ycoord \
--widthinpix --width=1 --oneelemstdout -q)
# Mask all the pixels that do not belong to the main object and set
# all the originally NaN-valued pixels to NaN again.
msk=$tmpdir/mask.fits
- astarithmetic $warpsnt $warplab $id ne nan where set-i \
+ astarithmetic $cropsnt $croplab $id ne nan where set-i \
i i $fillval eq nan where -g1 --output=$msk
# Set the masked image to the desired output of this step.
- mv $msk $warped_masked
+ mv $msk $cropped_masked
fi
# No '--segment' provided.
else
- warped_masked=$warped
+ cropped_masked=$cropped
+fi
+
+
+
+
+
+# Function to simplify checks in each dimension
+# ---------------------------------------------
+#
+# If the bottom-left edge of the crop is outside the image, the minimum
+# coordinate needs to be modified from the value in the 'ICF1PIX' keyword
+# of the output of the Crop program.
+first_pix_in_img () {
+
+ # This function takes a single argument: the dimension to work on.
+ dim=$1
+
+ # Find the minimim and maximum from the cropped image.
+ if [ $dim = x ]; then
+ width=$xwidthinpix
+ min=$(echo "$overlaprange" | awk '{print $1}')
+ max=$(echo "$overlaprange" | awk '{print $2}')
+ elif [ $dim = y ]; then
+ width=$ywidthinpix
+ min=$(echo "$overlaprange" | awk '{print $3}')
+ max=$(echo "$overlaprange" | awk '{print $4}')
+ else
+ cat <<EOF
+$scriptname: a bug! Please contact with us at bug-gnuastro@gnu.org. The value
of 'dim' is "$dim" and not recognized
+EOF
+ exit 1
+ fi
+
+ # Check if it is necessary to correct the position. If the bottom-left
+ # corner of the image (which defines the zero coordinate along both
+ # dimensions) has fallen outside of the image, then the minimum value
+ # of 'ICF1PIX' in that dimension will be '1'.
+ #
+ # Therefore, if the minimum value is not '1', then everything is fine
+ # and no correction is necessary! We'll just return the minimum
+ # value. Otherwise, we need to return a negative minimum pixel value to
+ # obtain the correct position of the object from 'ICF1PIX' (which only
+ # has the region of the input that overlapped with the output). It is
+ # necessary to add '1' because the first pixel has a coordiante of 1,
+ # not 0.
+ if [ $min = 1 ]; then
+ echo "$max $width" | awk '{if($1==$2) print 1; else print $1-$2+1}'
+ else
+ echo $min
+ fi
+}
+
+
+
+
+
+# In the FITS standard, the integer coordinates are defined on the center
+# of each pixel. On the other hand, the centers of objects can be anywhere
+# (not exactly in the center of the pixel!). In this scenario, we should
+# move the center of the object to the center of the pixel with a sub-pixel
+# warp. In this part of the script we will calculate the necessary
+# sub-pixel translation and use Warp to center the object on the pixel
+# grid. The user can disable this with the '--nocentering' option.
+if [ $nocentering = 0 ]; then
+
+ # Read the overlap range from the 'ICF1PIX' keyword (which is printed
+ # in all outputs of Crop).
+ overlaprange=$(astfits $cropped -h0 --keyvalue=ICF1PIX -q \
+ | sed -e's|:| |g' -e's|,| |')
+
+ # Calculate the position of the bottom-left pixel of the cropped image
+ # in relation to the input image.
+ minx=$(first_pix_in_img x)
+ miny=$(first_pix_in_img y)
+
+ # Due to cropping, the image coordinate should be shifted by a certain
+ # number (integer) of pixels (leaving the sub-pixel center intact).
+ stmcoord=$(echo "$xcenter $ycenter $minx $miny" \
+ | awk '{printf "%g %g\n", $1-$3+1, $2-$4+1}')
+
+ # Calculate the displacement (or value to give to '--translate' in
+ # Warp), and the new center of the star after translation (or value of
+ # '--center' in Crop):
+ # 1. At first we extract the sub-pixel displacement ('x' and 'y').
+ # 2. If the displacement is larger than 0.5, then after the warp, the
+ # center will be shifted by one pixel. Otherwise, if the
+ # displacement is smaller than 0.5, we shift the grid backwards by
+ # mutiplying it by -1 and don't add any shift.
+ # 3. Add 'xshift' and 'yshit' valut to integer of x and y coordinate
+ # and add final value to one.
+ warpcoord=$(echo "$stmcoord" \
+ | awk '{x=$1-int($1); y=$2-int($2); \
+ if(x>0.5){x=1-x; xshift=1} else {x*=-1; xshift=0}; \
+ if(y>0.5){y=1-y; yshift=1} else {y*=-1; yshift=0}; \
+ printf("%f,%f %d,%d\n", x, y, \
+ int($1)+xshift+1, int($2)+yshift+1)}')
+
+ # Warp image based on the measured displacement (first component of the
+ # output above).
+ warpped=$tmpdir/cropped-masked-warpforcenter.fits
+ DXY=$(echo "$warpcoord" | awk '{print $1}')
+ astwarp $cropped_masked --translate=$DXY --output=$warpped
+
+ # Crop image based on the calculated shift (second component of the
+ # output above).
+ centermsk=$tmpdir/cropped-masked-centered.fits
+ CXY=$(echo "$warpcoord" | awk '{print $2}')
+ astcrop $warpped -h1 \
+ --mode=img --output=$centermsk \
+ --center=$CXY --width=$xwidthinpix,$ywidthinpix
+else
+ # If the user did not want to correct the center of image, we'll use
+ # the raw crop as input for the next step.
+ centermsk=$cropped_masked
fi
@@ -786,7 +852,7 @@ if [ x"$normradiusmin" != x -a x"$normradiusmax" != x
]; then
if [ x"$sigmaclip" = x ]; then finalsigmaclip=""
else finalsigmaclip="--sigmaclip=$sigmaclip";
fi
- astscript-radial-profile $warped_masked --hdu=1 --rmax=$maxr \
+ astscript-radial-profile $centermsk --hdu=1 --rmax=$maxr \
--measure=$normop $finalsigmaclip \
--position-angle=$positionangle \
--tmpdir=$tmpdir --keeptmp \
@@ -827,7 +893,7 @@ fi
# but 1.868915 will be read as 'float64'. In the latter case, the output
# will become 'float64' also, which will cause problems later when we want
# to stack them together (some images will be 'float32', some 'float64').
-astarithmetic $warped_masked --hdu=1 $normvalue float32 / \
+astarithmetic $centermsk --hdu=1 $normvalue float32 / \
float32 --output=$output $quiet
diff --git a/bin/script/psf-subtract.sh b/bin/script/psf-subtract.sh
index 6f408240..395865ba 100755
--- a/bin/script/psf-subtract.sh
+++ b/bin/script/psf-subtract.sh
@@ -387,48 +387,23 @@ fi
-# Fake WCS info for no-WCS inputs and PSF image
-# ---------------------------------------------
-#
-# This fake WCS information is used if the input image does not have
-# astrometric information. In the same way, it is used for the PSF because
-# in principle the PSF image never has WCS information.
-
-# For the input image with no WCS info. The WCS injected is the following:
-# - crpix1,crpix2 is the center of the object to be PSF-modeled/subtracted.
-# - That pixel position is assigned to the crval1,crval2 sky position.
-# - The pixel scale is defined by cdelt1,cdelt2 (in arcsec).
-#
-# For PSF image, we assume it does not have WCS information. In this case:
-# - crpix1,crpix2 is the center of the image (i.e., the center of the PSF).
-# - That pixel position is assigned to the crval1,crval2 sky position.
-# - The pixel scale is defined by cdelt1,cdelt2 (in arcsec).
-
-crval1=180.0
-crval2=0.0
-cdelt1=1.0
-cdelt2=1.0
-cunits="deg,deg"
-
-
-
-
-
# Center coordinates and default object label
# -------------------------------------------
#
-# Obtain the coordinates of the center and generate a specific label for
-# the object consisting in its coordinates.
+# Obtain the coordinates of the center.
xcoord=$(echo "$center" | awk 'BEGIN{FS=","} {print $1}')
ycoord=$(echo "$center" | awk 'BEGIN{FS=","} {print $2}')
+
+# With the center coordinates, generate a specific label for the object
+# consisting in its coordinates.
objectid="$xcoord"_"$ycoord"
-# Define a temporal directory and the final output file
-# -----------------------------------------------------
+# Define a temporal directory and thefinal output file
+# ----------------------------------------------------
#
# Construct the temporary directory. If the user does not specify any
# directory, then a default one with the base name of the input image will
@@ -467,6 +442,43 @@ fi
+# 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)
+ xcenter=$(echo "$xycenter" | awk '{print $1}')
+ ycenter=$(echo "$xycenter" | awk '{print $2}')
+else
+ xcenter=$xcoord
+ ycenter=$ycoord
+fi
+
+
+
+
+
+# Scale the PSF with the given factor
+# -----------------------------------
+#
+# The PSF absolute pixel values are not relevant since they are used to be
+# normalized somehow. Here, the input PSF is scaled (multiplied) by the
+# specified factor (--fluxfactor) in order to appropiately obtain the PSF
+# brightness.
+psffluxscaled=$tmpdir/psf-flux-scaled-$objectid.fits
+astarithmetic $psf --hdu=$psfhdu $scale float32 x \
+ --output=$psffluxscaled $quiet
+
+
+
+
+
# Get the original inputs size
# ----------------------------
#
@@ -487,85 +499,17 @@ naxises=$(echo $axises \
| sed 's/n\/a//g' \
| awk '{print NF}')
if [ x$naxises = x4 ]; then
- naxis1=$(echo $axises | awk '{print $1}')
- naxix2=$(echo $axises | awk '{print $2}')
+ xaxis=$(echo $axises | awk '{print $1}')
+ yaxis=$(echo $axises | awk '{print $2}')
else
- naxis1=$(echo $axises | awk '{print $3}')
- naxis2=$(echo $axises | awk '{print $4}')
+ xaxis=$(echo $axises | awk '{print $3}')
+ yaxis=$(echo $axises | awk '{print $4}')
fi
-# 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
-
- # If the mode is WCS, the output should have WCS information
- wcsoutputinfo=""
- crval1=$xcoord
- crval2=$ycoord
- cdelt1=$(astfits $inputs --hdu $hdu --pixelscale --quiet | awk '{print
$1}')
- cdelt2=$(astfits $inputs --hdu $hdu --pixelscale --quiet | awk '{print
$2}')
-
-
-else
- # Fake WCS for the image is constructed by setting the center (in IMG)
- # to the center of the WCS information.
- center1=$xcoord
- center2=$ycoord
-
- fake_wcs=$tmpdir/fake-wcs-$objectid.fits
- echo "1 $center1 $center2 1 1 1 1 1 1 1" \
- | astmkprof --oversample=1 \
- --type=uint8 \
- --cunit=$cunits \
- --output=$fake_wcs \
- --cdelt=$cdelt1,$cdelt2 \
- --crval=$crval1,$crval2 \
- --crpix=$center1,$center2 \
- --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}')
-
- # If the mode is IMG, the output should not have WCS information
- wcsoutputinfo="--wcsfile=none"
-
-fi
-
-
-
-
# Get the center of the PSF
# -------------------------
#
@@ -579,68 +523,71 @@ psfnaxises=$(echo $psfaxises \
| sed 's/n\/a//g' \
| awk '{print NF}')
if [ x$psfnaxises = x4 ]; then
- psfnaxis1=$(echo $psfaxises | awk '{print $1}')
- psfnaxis2=$(echo $psfaxises | awk '{print $2}')
+ xpsfaxis=$(echo $psfaxises | awk '{print $1}')
+ ypsfaxis=$(echo $psfaxises | awk '{print $2}')
else
- psfnaxis1=$(echo $psfaxises | awk '{print $3}')
- psfnaxis2=$(echo $psfaxises | awk '{print $4}')
+ xpsfaxis=$(echo $psfaxises | awk '{print $3}')
+ ypsfaxis=$(echo $psfaxises | awk '{print $4}')
fi
-
-
-
-
-# WCS for the PSF image
-# ---------------------
-#
-# We assume that 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 from the information
-# defined above. In short: Center is the center of the PSF, but assignetd
-# to the center of the object/star in WCS coordinates.
-psfcenter1=$(echo $psfnaxis1 | awk '{print $1/2+0.5}')
-psfcenter2=$(echo $psfnaxis2 | awk '{print $1/2+0.5}')
-
-# Create the fake WCS image. No data, just WCS.
-psf_fake_wcs=$tmpdir/psf-fake-wcs.fits
-echo "1 $psfcenter1 $psfcenter2 1 1 1 1 1 1 1" \
- | astmkprof --oversample=1 \
- --type=uint8 \
- --cunit=$cunits \
- --output=$psf_fake_wcs \
- --cdelt=$cdelt1,$cdelt2 \
- --crval=$crval1,$crval2 \
- --crpix=$psfcenter1,$psfcenter2 \
- --mergedsize=$psfnaxis1,$psfnaxis2
+# To obtain the center of the PSF (in pixels), just divide by 2 the size of
+# the PSF image.
+xpsfcenter=$(astarithmetic $xpsfaxis float32 2.0 / --quiet)
+ypsfcenter=$(astarithmetic $ypsfaxis float32 2.0 / --quiet)
-# Inject fake WCS and scale the PSF by the given factor
-# -----------------------------------------------------
+# Translate the PSF image
+# -----------------------
#
-# The PSF absolute pixel values are not relevant since they are used to be
-# normalized somehow. Here, the input PSF is scaled (multiplied) by the
-# specified factor (--fluxfactor) in order to appropiately obtain the PSF
-# brightness. The fake WCS information is also added into the header.
-psf_wcs_scaled=$tmpdir/psf-wcsed-scaled-$objectid.fits
-astarithmetic $psf --hdu=$psfhdu $scale float32 x \
- --wcsfile=$psf_fake_wcs --output=$psf_wcs_scaled
+# In order to allocate the PSF into the center coordinates provided by the
+# user, it is necessary to compute the appropiate offsets along the X and Y
+# axis. After that, the PSF image is warped using that offsets.
+xdiff=$(astarithmetic $xcenter float32 $xpsfcenter float32 - --quiet)
+ydiff=$(astarithmetic $ycenter float32 $ypsfcenter float32 - --quiet)
+psftranslated=$tmpdir/psf-translated-$objectid.fits
+astwarp $psffluxscaled --translate=$xdiff,$ydiff \
+ --output=$psftranslated $quiet
-# Warp the PSF to into the correct position
-# -----------------------------------------
+
+# Crop the PSF image
+# ------------------
#
-# Once the PSF contains the WCS information to be properly positioned and
-# it has been flux scaled, here it is warped. The output is the model of
-# the object in the right position and with the appropriate flux level.
-psf_warped=$tmpdir/psf-warped-$objectid.fits
-astwarp $psf_wcs_scaled --gridfile=$inputs_wcs \
- --gridhdu $hdu --output=$psf_warped $quiet
+# Once the PSF has been situated appropiately into the right position, it
+# is necessary to crop it with a sub-pixel precision as well as with the
+# same size than the original input image. Otherwise the subtraction of the
+# PSF or scattered light field model would not be possible. Here, this
+# cropping is done.
+xrange=$(echo "$xdiff $xaxis $xcenter" \
+ | awk '{if($1<0) \
+ {i=int($1); \
+ d=-1*i+1; \
+ min=d; max=$2+d-1;} \
+ else {min=1; max=$2} \
+ i=int($3); {min=min+1; max=max+1}
+ printf "%d:%d", min, max}')
+
+yrange=$(echo "$ydiff $yaxis $ycenter" \
+ | awk '{if($1<0) \
+ {i=int($1); \
+ d=-1*i+1; \
+ min=d; max=$2+d-1;} \
+ else {min=1; max=$2} \
+ i=int($3); {min=min+1; max=max+1}
+ printf "%d:%d", min, max}')
+
+# Once the necessary crop parameters have been computed, use the option
+# '--section' to get the PSF image with the correct size and the center
+# situated where the user has specified.
+psfcropped=$tmpdir/psf-cropped-$objectid.fits
+astcrop $psftranslated --mode=img --section=$xrange,$yrange \
+ --output=$psfcropped $quiet
@@ -660,14 +607,12 @@ astwarp $psf_wcs_scaled --gridfile=$inputs_wcs \
# stars (using the PSF) or the subtraction without having nan values in the
# final subtracted image.
if [ x"$modelonly" = x0 ]; then
- astarithmetic $psf_warped --hdu=1 set-i \
+ astarithmetic $psfcropped --hdu=1 set-i \
i i isblank 0 where set-psf \
- $inputs --hdu=$hdu psf - float32 \
- $wcsoutputinfo --output=$output $quiet
+ $inputs --hdu=$hdu psf - --output=$output $quiet
else
- astarithmetic $psf_warped --hdu=1 set-i \
- i i isblank 0 where float32 \
- $wcsoutputinfo --output=$output $quiet
+ astarithmetic $psfcropped --hdu=1 set-i \
+ i i isblank 0 where --output=$output $quiet
fi
@@ -683,4 +628,3 @@ fi
if [ $keeptmp = 0 ]; then
rm -r $tmpdir
fi
-
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 8746deb1..ea589c2b 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -6699,6 +6699,7 @@ 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.
By doing this, the script will generate the radial profile of the entire stamp.
In this particular step we set it to @code{--normradii=500,510}.
+We also use the @option{--nocentering} option to disable sub-pixel warping in
this phase (it is only relevant for the central part of the PSF).
Furthermore, since there are several stars, we iterate over each row of the
catalog using a while loop.
@example
@@ -6708,6 +6709,7 @@ $ asttable outer/67510-6-10.fits \
| while read -r ra dec mag; do
astscript-psf-stamp label/67510-seg.fits \
--mode=wcs \
+ --nocentering \
--center=$ra,$dec \
--normradii=500,510 \
--widthinpix=1000,1000 \
@@ -6750,6 +6752,7 @@ $ asttable outer/67510-6-10.fits \
| while read -r ra dec mag; do
astscript-psf-stamp label/67510-seg.fits \
--mode=wcs \
+ --nocentering \
--center=$ra,$dec \
--normradii=20,30 \
--widthinpix=1000,1000 \
@@ -6770,17 +6773,9 @@ Just recall that we need to specify the number of inputs
into the stacking opera
@example
$ numimgs=$(echo outer/stamps/*.fits | wc -w)
$ astarithmetic outer/stamps/*.fits $numimgs 3 0.2 sigclip-mean \
- -g1 --output=outer/stack.fits --wcsfile=none --writeall
+ -g1 --output=outer/stack.fits --wcsfile=none
@end example
-Now opent the output of the stack image with following command.
-
-@example
-$ astscript-fits-view outer/stack.fits
-@end example
-
-The first extension represents the stack's output, whereas the second
extension indicates the number of images included in the final stack in each
pixel due to the @option{--writeall} option.
-
@noindent
Did you notice the @option{--wcsfile=none} option above?
With it, the stacked image no longer has any WCS information.
@@ -6855,7 +6850,7 @@ $ asttable inner/67510-12-13.fits \
$ numimgs=$(echo inner/stamps/*.fits | wc -w)
$ astarithmetic inner/stamps/*.fits $numimgs 3 0.2 sigclip-mean \
- -g1 --output=inner/stack.fits --wcsfile=none --writeall
+ -g1 --output=inner/stack.fits --wcsfile=none
$ astscript-fits-view inner/stack.fits inner/stamps/*.fits
@end example
@@ -6983,22 +6978,16 @@ 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 compute the scaling and keep it into a plain text file.
-Also, remove the temporary directory used previously.
+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.
@example
-$ mkdir uniting
$ rm -r stack_psfmodelscalefactor
-$ astscript-psf-scale-factor outer/stack.fits \
- --psf=inner/stack.fits --center=501,501 --quiet \
- --mode=img --normradii=10,15 --output=uniting/scale.txt
-$ cat uniting/scale.txt
+$ scale=$(astscript-psf-scale-factor outer/stack.fits \
+ --psf=inner/stack.fits --center=501,501 \
+ --mode=img --normradii=10,15 --quiet)
+$ echo $scale
@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}).
@@ -7006,7 +6995,6 @@ The inner part is first scaled, and all the pixels of the
outer image within the
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):
@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
@@ -7030,15 +7018,13 @@ 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
-$ astscript-psf-scale-factor outer/stack.fits \
- --mode=img --normradii=50,60 \
- --output=uniting/scale-bad.txt --quiet \
- --psf=inner/stack.fits --center=501,501
+$ scale=$(astscript-psf-scale-factor outer/stack.fits \
+ --psf=inner/stack.fits --center=501,501 \
+ --mode=img --normradii=50,60 --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_bad --output=psf-bad.fits
+ --scale=$scale --output=psf-bad.fits
$ astscript-fits-view psf-bad.fits psf.fits --ds9scale=minmax \
--ds9extra="-scale limits 0 50 -match scale" \
@@ -7087,41 +7073,27 @@ 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) of the
fainter star that we have used to construct the outer part of the PSF (note the
@option{--tail 1}).
+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):
@example
$ mkdir single-star
-$ center=$(asttable outer/67510-6-10.fits --sort phot_g_mean_mag \
- --column=ra,dec --tail 1 \
- | awk '@{printf "%.8f,%.8f\n", $1, $2@}')
+$ center=$(asttable flat/67510-bright.fits --sort phot_g_mean_mag \
+ --column=ra,dec --head 1 \
+ | awk '@{printf "%s,%s", $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,
this parameter is saved into a plain text file:
+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:
@example
-$ 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@}')
+$ 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)
@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.
@@ -7137,157 +7109,6 @@ $ 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 part 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.
-The pixel scale of this image is 0.55 arcsec/pixel, so with the line below we
shift the image along the R.A. direction by +1 pixel.
-
-@example
-$ original_center=$center
-$ modified_center=$(echo $center \
- | awk -F',' '@{printf "%.8f,%.8f\n", \
- $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 @option{--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-shifts.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.
@@ -7309,9 +7130,8 @@ We had 8 images for the outer part, so the S/N has only
improved by a factor of
When we multiply the final stacked PSF with 19, we are also scaling up the
noise by that same factor (most importantly: in the outer most regions where
there is almost no signal).
So the stacked image's noise-level is @mymath{19/3=6.3} times larger than the
noise of the input image.
This terrible noise-level is what you clearly see as the footprint of the PSF.
-In this case, to address the issue, we require additional stars to create the
outer part of the PSF and resolve the problem.
-Let's use the commands below to subtract the faintest of the bright-stars
catalog (note the use of @option{--tail} when finding the central position).
+To confirm this, let's use the commands below to subtract the faintest of the
bright-stars catalog (note the use of @option{--tail} when finding the central
position).
You will notice that the scale factor (@mymath{\sim1.3}) is now smaller than 3.
So when we multiply the PSF with this factor, the PSF's noise level is lower
than our input image and we should not see any footprint like before.
Note also that we are using a larger zoom factor, because this star is smaller
in the image.
@@ -35589,6 +35409,9 @@ A certain width (specified by @option{--widthinpix} in
pixels).
Centered at the coordinate specified by the option @option{--center} (it can
be in image/pixel or WCS coordinates, see @option{--mode}).
If no center is specified, then it is assumed that the object of interest is
already in the center of the image.
@item
+If the given coordinate has sub-pixel elements (for example, pixel coordinates
1.234,4.567), the pixel grid of the output will be warped so your given
coordinate falls in the center of the central pixel of the final output.
+This is very important for building the central parts of the PSF, but not too
effective for the middle or outer parts (to speed up the program in such cases,
you can disable it with the @option{--nocentering} option).
+@item
Normalized ``normalized'' by the value computed within the ring around the
center (at a radial distance between the two radii specified by the option
@option{--normradii}).
If no normalization ring is considered, the output will not be normalized.
@end itemize
@@ -35624,10 +35447,21 @@ This option thus accepts only two values:
@option{img} or @option{wcs}.
@itemx --center=FLT,FLT
The central position of the object.
This option is used for placing the center of the stamp.
-This parameter is used in @ref{Warp} to center and crop the image.
+This parameter is used in @ref{Crop} to center and crop the image.
The positions along each dimension must be separated by a comma (@key{,}).
The units of the coordinates are read based on the value to the
@option{--mode} option, see the examples above.
+The given coordinate for the central value can have sub-pixel elements (for
example, it falls on coordinate 123.4,567.8 of the input image pixel grid).
+In such cases, after cropping, this script will use Gnuastro's @ref{Warp} to
shift (or translate) the pixel grid by @mymath{-0.4} pixels along the
horizontal and @mymath{1-0.8=0.2} pixels along the vertical.
+Finally the newly added pixels (due to the warping) will be trimmed to have
your desired coordinate exactly in the center of the central pixel of the
output.
+This is very important (critical!) when you are constructing the central part
of the PSF.
+But for the outer parts it is not too effective, so to avoid wasting time for
the warping, you can simply use @option{--nocentering} to disable it.
+
+@item -d
+@itemx --nocentering
+Do not do the sub-pixel centering to a new pixel grid.
+See the description of the @option{--center} option for more.
+
@item -W INT,INT
@itemx --widthinpix=INT,INT
Size (width) of the output image stamp in pixels.
@@ -35861,11 +35695,8 @@ $ asttable catalog.fits | while read -r ra dec mag; do
\
@end example
The input should be an image containing the star that you want to match in
flux with the PSF.
-The output will be a three numbers that will be printed on the command-line,
or saved into a file if specified with the option @option{--output}.
-These numbers are:
-1. the multiplicative factor to scale the PSF image (given to @option{--psf})
to match in flux with the given star (which is located in @option{--center}
coordinate of the input image).
-2. the re-calculated x-center position if the option @option{--xshifts} is
used.
-3. the re-calculated y-center position if the option @option{--xshifts} is
used.
+The output will be a single number that is printed on the command-line.
+That number is the multiplicative factor to scale the PSF image (given to
@option{--psf}) to match in flux with the given star (which is located in
@option{--center} coordinate of the input image).
The scale factor will be calculated within the ring of pixels specified by the
option @option{--normradii}.
All the pixels within this ring will be separated from both the PSF and input
images.
@@ -35908,6 +35739,18 @@ This parameter is passed to Gnuastro's Crop program
make a crop for further proc
The positions along each dimension must be separated by a comma (@key{,}).
The units of the coordinates are interpreted based on the value to the
@option{--mode} option (see below).
+The given coordinate for the central value can have sub-pixel elements (for
example, it falls on coordinate 123.4,567.8 of the input image pixel grid).
+In such cases, after cropping, this script will use Gnuastro's @ref{Warp} to
shift (or translate) the pixel grid by @mymath{-0.4} pixels along the
horizontal and @mymath{1-0.8=0.2} pixels along the vertical.
+Finally the newly added pixels (due to the warping) will be trimmed to have
your desired coordinate exactly in the center of the central pixel of the
output.
+This is very important (critical!) when you are constructing the central part
of the PSF.
+But for the very far outer parts it may not too effective (should be checked),
or the target object may have already been centered at the requested coordinate.
+In such cases, to avoid wasting time for the warping, you can simply use
@option{--nocentering} to disable sub-pixel centering.
+
+@item -d
+@itemx --nocentering
+Do not do the sub-pixel centering to a new pixel grid.
+See the description of the @option{--center} option for more.
+
@item -O STR
@itemx --mode=STR
Interpret the center position of the object (values given to
@option{--center}) in image or WCS coordinates.
@@ -35945,17 +35788,6 @@ The result is that all objects and clumps that
contaminate the central source ar
Sigma clipping parameters used in the end to find the final scale factor from
the distribution of all pixels used.
For more on sigma-clipping, see @ref{Sigma clipping}.
-@item --xshifts=FLT,FLT,FLT
-Shifts on the centering position provided by @option{--center} in units of
pixels.
-The three parameters used by the program @code{seq} to create the different
shifted positions.
-For example, @option{--xshifts=-3.0,0.4,+1.0} will shift the center position
from -3.0 pixels to +1.0 pixels in step of 0.4 pixels, along the x direction.
-It will compute the scaling factor for all of those position, and the scale
factor for wich the standard deviation is the smallest will be provided as the
best center.
-The re-computed center position is also provided as part of the output.
-
-@item --yshifts=FLT,FLT,FLT
-Same than the @option{--xshifts} but for the y direction.
-See above for more information.
-
@item -t
@itemx --tmpdir
Directory to keep temporary files during the execution of the script.
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master e985712d: psf-scale-factor, psf-stamp and psf-subtract: revered to eaab08,
Mohammad Akhlaghi <=